MolLib/Parameter.c

   1 /*
   2  *  Parameter.c
   3  *
   4  *  Created by Toshi Nagata on 06/03/11.
   5  *  Copyright 2006-2008 Toshi Nagata. All rights reserved.
   6  *
   7  This program is free software; you can redistribute it and/or modify
   8  it under the terms of the GNU General Public License as published by
   9  the Free Software Foundation version 2 of the License.
  10
  11  This program is distributed in the hope that it will be useful,
  12  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14  GNU General Public License for more details.
  15  */
  16
  17 #include "MolLib.h"
  18 #include <string.h>
  19 #include <ctype.h>
  20 #include <stdarg.h>
  21 #include <math.h>
  22
  23 /*  Global parameter: it is initialized by the first call to ParameterReadFromFile()  */
  24 Parameter *gBuiltinParameters = NULL;
  25
  26 /*  Global parameter  */
  27 ElementPar *gElementParameters = NULL;
  28 Int gCountElementParameters = 0;
  29
  30 static Parameter *sParameterRoot = NULL;
  31 static int sParameterUntitledCount = 0;
  32
  33 /*  Global struct for storing information for global parameter files  */
  34 GlobalParInfoRecord gGlobalParInfo;
  35
  36 #pragma mark ====== Parameter Alloc/Release ======
  37
  38 Parameter *
  39 ParameterNew(void)
  40 {
  41         char name[40];
  42         Parameter *par = (Parameter *)calloc(sizeof(Parameter), 1);
  43         if (par == NULL)
  44                 Panic("Cannot allocate new parameter record");
  45         snprintf(name, sizeof name, "Untitled %d", sParameterUntitledCount++);
  46         ObjectInit((Object *)par, (Object **)&sParameterRoot, name);
  47         return par;
  48 }
  49
  50 Parameter *
  51 ParameterDuplicate(const Parameter *par)
  52 {
  53         Parameter *npar = ParameterNew();
  54         if (par->bondPars != NULL) {
  55                 npar->bondPars = (BondPar *)malloc(sizeof(BondPar) * par->nbondPars);
  56                 if (npar->bondPars == NULL)
  57                         goto release;
  58                 memmove(npar->bondPars, par->bondPars, sizeof(BondPar) * par->nbondPars);
  59                 npar->nbondPars = par->nbondPars;
  60         }
  61         if (par->anglePars != NULL) {
  62                 npar->anglePars = (AnglePar *)malloc(sizeof(AnglePar) * par->nanglePars);
  63                 if (npar->anglePars == NULL)
  64                         goto release;
  65                 memmove(npar->anglePars, par->anglePars, sizeof(AnglePar) * par->nanglePars);
  66                 npar->nanglePars = par->nanglePars;
  67         }
  68         if (par->dihedralPars != NULL) {
  69                 npar->dihedralPars = (TorsionPar *)malloc(sizeof(TorsionPar) * par->ndihedralPars);
  70                 if (npar->dihedralPars == NULL)
  71                         goto release;
  72                 memmove(npar->dihedralPars, par->dihedralPars, sizeof(TorsionPar) * par->ndihedralPars);
  73                 npar->ndihedralPars = par->ndihedralPars;
  74         }
  75         if (par->improperPars != NULL) {
  76                 npar->improperPars = (TorsionPar *)malloc(sizeof(TorsionPar) * par->nimproperPars);
  77                 if (npar->improperPars == NULL)
  78                         goto release;
  79                 memmove(npar->improperPars, par->improperPars, sizeof(TorsionPar) * par->nimproperPars);
  80                 npar->nimproperPars = par->nimproperPars;
  81         }
  82         if (par->vdwPars != NULL) {
  83                 npar->vdwPars = (VdwPar *)malloc(sizeof(VdwPar) * par->nvdwPars);
  84                 if (npar->vdwPars == NULL)
  85                         goto release;
  86                 memmove(npar->vdwPars, par->vdwPars, sizeof(VdwPar) * par->nvdwPars);
  87                 npar->nvdwPars = par->nvdwPars;
  88         }
  89         if (par->vdwpPars != NULL) {
  90                 npar->vdwpPars = (VdwPairPar *)malloc(sizeof(VdwPairPar) * par->nvdwpPars);
  91                 if (npar->vdwpPars == NULL)
  92                         goto release;
  93                 memmove(npar->vdwpPars, par->vdwpPars, sizeof(VdwPairPar) * par->nvdwpPars);
  94                 npar->nvdwpPars = par->nvdwpPars;
  95         }
  96         if (par->vdwCutoffPars != NULL) {
  97                 npar->vdwCutoffPars = (VdwCutoffPar *)malloc(sizeof(VdwCutoffPar) * par->nvdwCutoffPars);
  98                 if (npar->vdwCutoffPars == NULL)
  99                         goto release;
 100                 memmove(npar->vdwCutoffPars, par->vdwCutoffPars, sizeof(VdwCutoffPar) * par->nvdwCutoffPars);
 101                 npar->nvdwCutoffPars = par->nvdwCutoffPars;
 102         }
 103 /*      if (par->atomPars != NULL) {
 104                 npar->atomPars = (ElementPar *)malloc(sizeof(ElementPar) * par->natomPars);
 105                 if (npar->atomPars == NULL)
 106                         goto release;
 107                 memmove(npar->atomPars, par->atomPars, sizeof(ElementPar) * par->natomPars);
 108                 npar->natomPars = par->natomPars;
 109         } */
 110         return npar;
 111 release:
 112         ParameterRelease(npar);
 113         return NULL;
 114 }
 115
 116 Parameter *
 117 ParameterWithName(const char *name)
 118 {
 119         return (Parameter *)ObjectWithName(name, (Object *)sParameterRoot);
 120 }
 121
 122 /*  Assign a unique name to this parameter record  */
 123 void
 124 ParameterSetName(Parameter *par, const char *name)
 125 {
 126         ObjectSetName((Object *)par, name, (Object *)sParameterRoot);
 127 }
 128
 129 const char *
 130 ParameterGetName(Parameter *par)
 131 {
 132         return ObjectGetName((Object *)par);
 133 }
 134
 135 void
 136 ParameterRetain(Parameter *par)
 137 {
 138         ObjectIncrRefCount((Object *)par);
 139 }
 140
 141 void
 142 ParameterRelease(Parameter *par)
 143 {
 144         if (ObjectDecrRefCount((Object *)par) == 0) {
 145                 if (par->bondPars != NULL)
 146                         free(par->bondPars);
 147                 if (par->anglePars != NULL)
 148                         free(par->anglePars);
 149                 if (par->dihedralPars != NULL)
 150                         free(par->dihedralPars);
 151                 if (par->improperPars != NULL)
 152                         free(par->improperPars);
 153                 if (par->vdwPars != NULL)
 154                         free(par->vdwPars);
 155                 if (par->vdwpPars != NULL)
 156                         free(par->vdwpPars);
 157                 if (par->vdwCutoffPars != NULL)
 158                         free(par->vdwCutoffPars);
 159         /*      if (par->atomPars != NULL)
 160                         free(par->atomPars); */
 161                 ObjectDealloc((Object *)par, (Object **)&sParameterRoot);
 162         }
 163 }
 164
 165 #pragma mark ====== ParameterRef Definitions ======
 166
 167 ParameterRef *
 168 ParameterRefNew(struct Molecule *mol, int type, int idx)
 169 {
 170         ParameterRef *pref = (ParameterRef *)calloc(sizeof(ParameterRef), 1);
 171         if (pref != NULL) {
 172                 pref->mol = mol;
 173                 if (mol != NULL)
 174                         MoleculeRetain(mol);
 175                 pref->parType = type;
 176                 pref->idx = idx;
 177         }
 178         return pref;
 179 }
 180
 181 void
 182 ParameterRefRelease(ParameterRef *pref)
 183 {
 184         if (pref != NULL) {
 185                 if (pref->mol != NULL)
 186                         MoleculeRelease(pref->mol);
 187                 free(pref);
 188         }
 189 }
 190
 191 UnionPar *
 192 ParameterGetUnionParFromTypeAndIndex(Parameter *par, int type, int index)
 193 {
 194         if (par == NULL)
 195                 return NULL;
 196         switch (type) {
 197                 case kBondParType:
 198                         if (index >= 0 && index < par->nbondPars)
 199                                 return (UnionPar *)(par->bondPars + index);
 200                         else return NULL;
 201                 case kAngleParType:
 202                         if (index >= 0 && index < par->nanglePars)
 203                                 return (UnionPar *)(par->anglePars + index);
 204                         else return NULL;
 205                 case kDihedralParType:
 206                         if (index >= 0 && index < par->ndihedralPars)
 207                                 return (UnionPar *)(par->dihedralPars + index);
 208                         else return NULL;
 209                 case kImproperParType:
 210                         if (index >= 0 && index < par->nimproperPars)
 211                                 return (UnionPar *)(par->improperPars + index);
 212                         else return NULL;
 213                 case kVdwParType:
 214                         if (index >= 0 && index < par->nvdwPars)
 215                                 return (UnionPar *)(par->vdwPars + index);
 216                         else return NULL;
 217                 case kVdwPairParType:
 218                         if (index >= 0 && index < par->nvdwpPars)
 219                                 return (UnionPar *)(par->vdwpPars + index);
 220                         else return NULL;
 221                 case kVdwCutoffParType:
 222                         if (index >= 0 && index < par->nvdwCutoffPars)
 223                                 return (UnionPar *)(par->vdwCutoffPars + index);
 224                         else return NULL;
 225         /*      case kElementParType:
 226                         if (index >= 0 && index < par->natomPars)
 227                                 return (UnionPar *)(par->atomPars + index);
 228                         else return NULL; */
 229                 default:
 230                         return NULL;
 231         }
 232 }
 233
 234 Int
 235 ParameterGetCountForType(Parameter *par, int type)
 236 {
 237         if (par == NULL)
 238                 return 0;
 239         switch (type) {
 240                 case kBondParType:
 241                         return par->nbondPars;
 242                 case kAngleParType:
 243                         return par->nanglePars;
 244                 case kDihedralParType:
 245                         return par->ndihedralPars;
 246                 case kImproperParType:
 247                         return par->nimproperPars;
 248                 case kVdwParType:
 249                         return par->nvdwPars;
 250                 case kVdwPairParType:
 251                         return par->nvdwpPars;
 252                 case kVdwCutoffParType:
 253                         return par->nvdwCutoffPars;
 254         /*      case kElementParType:
 255                         return par->natomPars; */
 256                 default:
 257                         return 0;
 258         }
 259 }
 260
 261 Int
 262 ParameterGetSizeForType(int type)
 263 {
 264         switch (type) {
 265                 case kBondParType:
 266                         return sizeof(BondPar);
 267                 case kAngleParType:
 268                         return sizeof(AnglePar);
 269                 case kDihedralParType:
 270                         return sizeof(TorsionPar);
 271                 case kImproperParType:
 272                         return sizeof(TorsionPar);
 273                 case kVdwParType:
 274                         return sizeof(VdwPar);
 275                 case kVdwPairParType:
 276                         return sizeof(VdwPairPar);
 277                 case kVdwCutoffParType:
 278                         return sizeof(VdwCutoffPar);
 279                         /*      case kElementParType:
 280                          return par->natomPars; */
 281                 default:
 282                         return 0;
 283         }
 284 }
 285
 286 void
 287 ParameterCopyOneWithType(UnionPar *dst, const UnionPar *src, int type)
 288 {
 289         int size = ParameterGetSizeForType(type);
 290         memmove(dst, src, size);
 291 }
 292
 293 UnionPar *
 294 ParameterRefGetPar(ParameterRef *pref)
 295 {
 296         Parameter *par;
 297         if (pref == NULL)
 298                 return NULL;
 299         if (pref->mol != NULL)
 300                 par = pref->mol->par;
 301         else par = gBuiltinParameters;
 302         if (par == NULL)
 303                 return NULL;
 304         return ParameterGetUnionParFromTypeAndIndex(par, pref->parType, pref->idx);
 305 }
 306
 307 #pragma mark ====== Insert/Delete (for MolAction) ======
 308
 309 int
 310 ParameterInsert(Parameter *par, Int type, const UnionPar *up, struct IntGroup *where)
 311 {
 312         Int i, n1, n2, size, *ip;
 313         void *p;
 314         if (par == NULL)
 315                 return -1;
 316         switch (type) {
 317                 case kBondParType:
 318                         p = &par->bondPars;
 319                         ip = &par->nbondPars;
 320                         size = sizeof(BondPar);
 321                         break;
 322                 case kAngleParType:
 323                         p = &par->anglePars;
 324                         ip = &par->nanglePars;
 325                         size = sizeof(AnglePar);
 326                         break;
 327                 case kDihedralParType:
 328                         p = &par->dihedralPars;
 329                         ip = &par->ndihedralPars;
 330                         size = sizeof(TorsionPar);
 331                         break;
 332                 case kImproperParType:
 333                         p = &par->improperPars;
 334                         ip = &par->nimproperPars;
 335                         size = sizeof(TorsionPar);
 336                         break;
 337                 case kVdwParType:
 338                         p = &par->vdwPars;
 339                         ip = &par->nvdwPars;
 340                         size = sizeof(VdwPar);
 341                         break;
 342                 case kVdwPairParType:
 343                         p = &par->vdwpPars;
 344                         ip = &par->nvdwpPars;
 345                         size = sizeof(VdwPairPar);
 346                         break;
 347                 case kVdwCutoffParType:
 348                         p = &par->vdwCutoffPars;
 349                         ip = &par->nvdwCutoffPars;
 350                         size = sizeof(VdwCutoffPar);
 351                         break;
 352         /*      case kElementParType:
 353                         p = &par->atomPars;
 354                         ip = &par->natomPars;
 355                         size = sizeof(ElementPar);
 356                         break; */
 357                 default:
 358                         return -1;
 359         }
 360         n2 = 0;
 361         for (i = 0; (n1 = IntGroupGetNthPoint(where, i)) >= 0; i++) {
 362                 if (InsertArray(p, ip, size, n1, 1, (up ? up + i : NULL)) == NULL)
 363                         return n2;
 364                 n2++;
 365         }
 366         return n2;
 367 }
 368
 369 static int
 370 sParameterDeleteOrCopy(Parameter *par, Int type, UnionPar *up, struct IntGroup *where, Int copyflag)
 371 {
 372         Int i, n1, n2, size, *ip;
 373         char **p;
 374         if (par == NULL)
 375                 return -1;
 376         switch (type) {
 377                 case kBondParType:
 378                         p = (char **)&par->bondPars;
 379                         ip = &par->nbondPars;
 380                         size = sizeof(BondPar);
 381                         break;
 382                 case kAngleParType:
 383                         p = (char **)&par->anglePars;
 384                         ip = &par->nanglePars;
 385                         size = sizeof(AnglePar);
 386                         break;
 387                 case kDihedralParType:
 388                         p = (char **)&par->dihedralPars;
 389                         ip = &par->ndihedralPars;
 390                         size = sizeof(TorsionPar);
 391                         break;
 392                 case kImproperParType:
 393                         p = (char **)&par->improperPars;
 394                         ip = &par->nimproperPars;
 395                         size = sizeof(TorsionPar);
 396                         break;
 397                 case kVdwParType:
 398                         p = (char **)&par->vdwPars;
 399                         ip = &par->nvdwPars;
 400                         size = sizeof(VdwPar);
 401                         break;
 402                 case kVdwPairParType:
 403                         p = (char **)&par->vdwpPars;
 404                         ip = &par->nvdwpPars;
 405                         size = sizeof(VdwPairPar);
 406                         break;
 407                 case kVdwCutoffParType:
 408                         p = (char **)&par->vdwCutoffPars;
 409                         ip = &par->nvdwCutoffPars;
 410                         size = sizeof(VdwCutoffPar);
 411                         break;
 412         /*      case kElementParType:
 413                         p = (char **)&par->atomPars;
 414                         ip = &par->natomPars;
 415                         size = sizeof(ElementPar);
 416                         break; */
 417                 default:
 418                         return -1;
 419         }
 420         n2 = 0;
 421         for (i = IntGroupGetCount(where) - 1; i >= 0 && (n1 = IntGroupGetNthPoint(where, i)) >= 0; i--) {
 422                 if (n1 >= *ip)
 423                         return -1; /*  Internal error  */
 424                 if (copyflag) {
 425                         if (up != NULL)
 426                                 memmove(up + i, *p + size * n1, size);
 427                 } else {
 428                         if (DeleteArray(p, ip, size, n1, 1, (up ? up + i : NULL)) == NULL)
 429                                 return n2;
 430                 }
 431                 n2++;
 432         }
 433         return n2;
 434 }
 435
 436 int
 437 ParameterDelete(Parameter *par, Int type, UnionPar *up, struct IntGroup *where)
 438 {
 439         return sParameterDeleteOrCopy(par, type, up, where, 0);
 440 }
 441
 442 int
 443 ParameterCopy(Parameter *par, Int type, UnionPar *up, struct IntGroup *where)
 444 {
 445         return sParameterDeleteOrCopy(par, type, up, where, 1);
 446 }
 447
 448 /*  Renumber the atom index field according to the conversion table.  If the atom index
 449     points to a non-existing atom, returns non-zero.  */
 450 int
 451 ParameterRenumberAtoms(Int type, UnionPar *up, Int oldnatoms, const Int *old2new)
 452 {
 453 #define RENUMBER_FIELD(_tp, _field) \
 454   (((_tp *)up)->_field >= 0 && ((_tp *)up)->_field < oldnatoms && \
 455         (old2new[((_tp *)up)->_field] >= 0 ? \
 456       ((((_tp *)up)->_field = old2new[((_tp *)up)->_field]), 0) : \
 457       ((((_tp *)up)->_field = kAtomTypeWildcard), 1)))
 458         switch (type) {
 459                 case kBondParType:
 460                         return RENUMBER_FIELD(BondPar, type1) + RENUMBER_FIELD(BondPar, type2);
 461                 case kAngleParType:
 462                         return RENUMBER_FIELD(AnglePar, type1) + RENUMBER_FIELD(AnglePar, type2) + RENUMBER_FIELD(AnglePar, type3);
 463                 case kDihedralParType:
 464                 case kImproperParType:
 465                         return RENUMBER_FIELD(TorsionPar, type1) + RENUMBER_FIELD(TorsionPar, type2) + RENUMBER_FIELD(TorsionPar, type3) + RENUMBER_FIELD(TorsionPar, type4);
 466                 case kVdwParType:
 467                         return RENUMBER_FIELD(VdwPar, type1);
 468                 case kVdwPairParType:
 469                         return RENUMBER_FIELD(VdwPairPar, type1) + RENUMBER_FIELD(VdwPairPar, type2);
 470                 case kVdwCutoffParType:
 471                         if (((VdwCutoffPar *)up)->type == 1)
 472                                 return RENUMBER_FIELD(VdwCutoffPar, n1) + RENUMBER_FIELD(VdwCutoffPar, n2) + RENUMBER_FIELD(VdwCutoffPar, n3) || RENUMBER_FIELD(VdwCutoffPar, n4);
 473                         else return 0;
 474                 default:
 475                         return -1;
 476         }
 477 }
 478
 479 #pragma mark ====== Load from Files ======
 480
 481 static int
 482 s_AppendWarning(char **ptr, const char *fmt, ...)
 483 {
 484         char *buf;
 485         int len, n, nn;
 486         va_list ap;
 487         if (ptr == NULL)
 488                 return 0;
 489         if (*ptr == NULL) {
 490                 *ptr = (char *)malloc(128);
 491                 if (*ptr == NULL)
 492                         return -1;
 493                 (*ptr)[0] = 0;
 494                 len = 128;
 495                 n = 0;
 496         } else {
 497                 n = strlen(*ptr);
 498                 for (len = 128; len <= n; len *= 2);
 499         }
 500         va_start(ap, fmt);
 501         nn = vasprintf(&buf, fmt, ap);
 502         if (nn < 0)
 503                 return nn;
 504         if (len <= n + nn) {
 505                 while (len <= n + nn)
 506                         len *= 2;
 507                 *ptr = (char *)realloc(*ptr, len);
 508                 if (*ptr == NULL)
 509                         return -1;
 510         }
 511         strncpy(*ptr + n, buf, nn);
 512         *(*ptr + n + nn) = 0;
 513         free(buf);
 514         return nn;
 515 }
 516
 517 int
 518 ParameterGlobalParIndexForSrcIndex(int src)
 519 {
 520         int i;
 521         if (src == 0 || src == -1)
 522                 return src - 1;
 523         for (i = 0; i < gGlobalParInfo.count; i++) {
 524                 if (gGlobalParInfo.files[i].src == src)
 525                         return i;
 526         }
 527         return -2;
 528 }
 529
 530 int
 531 ParameterCommentIndexForGlobalFileName(const char *p)
 532 {
 533         const char *pp = p + strlen(p), *p1 = NULL, *p2 = NULL;
 534         char buf[128];
 535         while (--pp >= p) {
 536                 if (p2 == NULL && *pp == '.')
 537                         p2 = pp;
 538                 if (p1 == NULL && (*pp == '\'' || *pp == '/'))
 539                         p1 = pp + 1;
 540         }
 541         if (p1 == NULL)
 542                 p1 = p;
 543         if (p2 == NULL)
 544                 p2 = p + strlen(p);
 545         snprintf(buf, sizeof(buf), "%.*s", (int)(p2 - p1), p1);
 546         return ParameterCommentIndex(buf);
 547 }
 548
 549 int
 550 ParameterCompare(const UnionPar *up1, const UnionPar *up2, int type)
 551 {
 552         switch (type) {
 553                 case kBondParType: {
 554                         const BondPar *bp1 = &up1->bond;
 555                         const BondPar *bp2 = &up2->bond;
 556                         return (((bp1->type1 == bp2->type1 && bp1->type2 == bp2->type2)
 557                                          || (bp1->type1 == bp2->type2 && bp1->type2 == bp2->type1))
 558                                         && fabs(bp1->k - bp2->k) < 1e-8 && fabs(bp1->r0 - bp2->r0) < 1e-8);
 559                 }
 560                 case kAngleParType: {
 561                         const AnglePar *ap1 = &up1->angle;
 562                         const AnglePar *ap2 = &up2->angle;
 563                         return (ap1->type2 == ap2->type2
 564                                         && ((ap1->type1 == ap2->type1 && ap1->type3 == ap2->type3)
 565                                                 || (ap1->type1 == ap2->type3 && ap1->type3 == ap2->type1))
 566                                         && fabs(ap1->k - ap2->k) < 1e-8 && fabs(ap1->a0 - ap2->a0) < 1e-8);
 567                 }
 568                 case kDihedralParType:
 569                 case kImproperParType: {
 570                         const TorsionPar *tp1 = &up1->torsion;
 571                         const TorsionPar *tp2 = &up2->torsion;
 572                         int i;
 573                         if (tp1->mult != tp2->mult)
 574                                 return 0;
 575                         if (type == kDihedralParType) {
 576                                 if ((tp1->type1 != tp2->type1 || tp1->type2 != tp2->type2 || tp1->type3 != tp2->type3 || tp1->type4 != tp2->type4)
 577                                         && (tp1->type1 != tp2->type4 || tp1->type2 != tp2->type3 || tp1->type3 != tp2->type2 || tp1->type4 != tp2->type1))
 578                                         return 0;
 579                         } else {
 580                                 if (tp1->type3 != tp2->type3)
 581                                         return 0;
 582                                 if ((tp1->type1 != tp2->type1 || tp1->type2 != tp2->type2 || tp1->type4 != tp2->type4)
 583                                         && (tp1->type1 != tp2->type1 || tp1->type2 != tp2->type4 || tp1->type4 != tp2->type2)
 584                                         && (tp1->type1 != tp2->type2 || tp1->type2 != tp2->type1 || tp1->type4 != tp2->type4)
 585                                         && (tp1->type1 != tp2->type2 || tp1->type2 != tp2->type4 || tp1->type4 != tp2->type1)
 586                                         && (tp1->type1 != tp2->type4 || tp1->type2 != tp2->type1 || tp1->type4 != tp2->type2)
 587                                         && (tp1->type1 != tp2->type4 || tp1->type2 != tp2->type2 || tp1->type4 != tp2->type1))
 588                                         return 0;
 589                         }
 590                         for (i = 0; i < tp1->mult; i++) {
 591                                 if (tp1->period[i] != tp2->period[i] || fabs(tp1->k[i] - tp2->k[i]) > 1e-8 || fabs(tp1->phi0[i] - tp2->phi0[i]) > 1e-8)
 592                                         return 0;
 593                         }
 594                         return 1;
 595                 }
 596                 case kVdwParType: {
 597                         const VdwPar *vp1 = &up1->vdw;
 598                         const VdwPar *vp2 = &up2->vdw;
 599                         return (vp1->type1 == vp2->type1 && fabs(vp1->eps - vp2->eps) < 1e-6 && fabs(vp1->r_eq - vp2->r_eq) < 1e-6 && fabs(vp1->eps14 - vp2->eps14) < 1e-4 && fabs(vp1->r_eq14 - vp2->r_eq14) < 1e-4);
 600                 }
 601                 case kVdwPairParType: {
 602                         const VdwPairPar *vp1 = &up1->vdwp;
 603                         const VdwPairPar *vp2 = &up2->vdwp;
 604                         return (vp1->type1 == vp2->type1 && fabs(vp1->eps - vp2->eps) < 1e-6 && fabs(vp1->r_eq - vp2->r_eq) < 1e-6 && fabs(vp1->eps14 - vp2->eps14) < 1e-4 && fabs(vp1->r_eq14 - vp2->r_eq14) < 1e-4);
 605                 }
 606                 case kVdwCutoffParType: {
 607                         const VdwCutoffPar *vp1 = &up1->vdwcutoff;
 608                         const VdwCutoffPar *vp2 = &up2->vdwcutoff;
 609                         if (vp1->type != vp2->type)
 610                                 return 0;
 611                         if (vp1->type == 0) {
 612                                 if ((vp1->n1 != vp2->n1 || vp1->n2 != vp2->n2) && (vp1->n1 != vp2->n2 || vp1->n2 != vp2->n1))
 613                                         return 0;
 614                         } else {
 615                                 if (vp1->n1 != vp2->n1 || vp1->n2 != vp2->n2 || vp1->n3 != vp2->n3 || vp1->n4 != vp2->n4)
 616                                         return 0;
 617                         }
 618                         return fabs(vp1->cutoff - vp2->cutoff) < 1e-8;
 619                 }
 620         }
 621         return 0;
 622 }
 623
 624 /*  bp can also be AnglePar *, TorsionPar *, etc.  */
 625 static void
 626 s_StoreComment(int parType, BondPar *bp, char *p, const char *src)
 627 {
 628         char *s, *pp, *pcom, buf[40];
 629         int embedded = 0, src_idx;
 630         if (p == NULL)
 631                 return;
 632         while (isspace(*p) || *p == ';' || *p == '!')
 633                 p++;
 634         pcom = p;
 635         if (src == NULL && *p == '[') {
 636                 /*  Source is embedded  */
 637                 int len;
 638                 s = p + 1;
 639                 p += 2;
 640                 while (*p != ']' && *p != 0)
 641                         p++;
 642                 len = p - s;
 643                 if (len >= sizeof(buf))
 644                         len = sizeof(buf) - 1;
 645                 strncpy(buf, s, len);
 646                 buf[len] = 0;
 647                 src = buf;
 648                 if (*p != 0) {
 649                         p++;
 650                         while (isspace(*p))
 651                                 p++;
 652                 }
 653                 embedded = 1;
 654         }
 655         pp = p;
 656         while (*pp != 0 && *pp != '\r' && *pp != '\n')
 657                 pp++;
 658         *pp = 0;
 659         if (src != NULL && *src != 0) {
 660                 src_idx = ParameterCommentIndex(src);
 661                 if (embedded) {
 662                         /*  Compare with already known global parameters, and src is set if the same one is found  */
 663                         int i;
 664                         UnionPar *up2;
 665                         for (i = 0; (up2 = ParameterGetUnionParFromTypeAndIndex(gBuiltinParameters, parType, i)) != NULL; i++) {
 666                                 if (up2->bond.src == src_idx && ParameterCompare((UnionPar *)bp, up2, parType)) {
 667                                         bp->src = src_idx;
 668                                         break;
 669                                 }
 670                         }
 671                         if (up2 == NULL) {
 672                                 /*  Not found: embedded source is retained, and this entry is regarded as "no source"  */
 673                                 bp->src = 0;
 674                                 p = pcom;
 675                         }
 676                 } else bp->src = src_idx;
 677         }
 678         if (*p != 0)
 679                 bp->com = ParameterCommentIndex(p);
 680 }
 681
 682 /*  bp can also be BondPar*, AnglePar *, TorsionPar *, etc.  */
 683 static void
 684 s_CommentToString(char *buf, int bufsize, void *bp)
 685 {
 686         const char *src, *com;
 687         src = ParameterGetComment(((BondPar *)bp)->src);
 688         com = ParameterGetComment(((BondPar *)bp)->com);
 689         if (src == NULL && com == NULL)
 690                 buf[0] = 0;
 691         else if (src != NULL)
 692                 snprintf(buf, bufsize, " ![%s] %s", src, (com == NULL ? "" : com));
 693         else
 694                 snprintf(buf, bufsize, " ! %s", com);
 695 }
 696
 697 int
 698 ParameterReadFromString(Parameter *par, char *buf, char **wbufp, const char *fname, int lineNumber, int src_idx)
 699 {
 700         int i, len, options;
 701         char com[12], com2[12], type[4][8];
 702         Int itype[4];
 703         float val[6];  /*  not Double  */
 704         Int ival[12];
 705         int n;
 706         int retval = 0;
 707         const char *src;
 708         if (sscanf(buf, " %11s", com) <= 0 || !isalpha(com[0]))
 709                 return 0;
 710         len = strlen(com);
 711         for (i = 0; i < len; i++)
 712                 com[i] = tolower(com[i]);
 713         if (strncmp(com, "include", len) == 0) {
 714                 char c, *p, *pp;
 715                 int wc1;
 716                 char *wbuf1;
 717                 i = 0;
 718                 for ( ; (c = buf[len]) != 0; len++) {
 719                         if (c == ' ' || c == '\t')
 720                                 continue;
 721                         if (c == '\"') {
 722                                 if (i == 0)
 723                                         continue;
 724                                 else break;
 725                         }
 726                         if (c == '\\') {
 727                                 len++;
 728                                 buf[i] = buf[len];
 729                                 if (buf[len] == 0)
 730                                         break;
 731                                 i++;
 732                                 continue;
 733                         }
 734                         buf[i++] = c;
 735                 }
 736                 buf[i] = 0;
 737                 p = (char *)malloc(strlen(fname) + i + 1);
 738                 if (p == NULL) {
 739                         return -1;
 740                 }
 741                 strcpy(p, fname);
 742                 pp = strrchr(p, '/');
 743                 if (pp == NULL)
 744                         strcpy(p, buf);
 745                 else
 746                         strcpy(p + 1, buf);
 747                 i = ParameterReadFromFile(par, p, &wbuf1, &wc1);
 748                 if (wbuf1 != NULL) {
 749                         s_AppendWarning(wbufp, "In included file %s:\n%s", p, wbuf1);
 750                         free(wbuf1);
 751                         retval = wc1;
 752                 }
 753                 free(p);
 754                 if (i < 0) {
 755                         retval = i;
 756                 }
 757                 return retval;
 758         }
 759         if (par == gBuiltinParameters && src_idx == 0) {
 760                 /*  Comes here only when the reading file is "default.par" at the initialization of the built-in parameters. */
 761                 /*  In this case, only "include" statements are allowed.  */
 762                 return -1;
 763         }
 764         if (src_idx == 0)
 765                 src = NULL;
 766         else src = ParameterGetComment(src_idx);
 767         options = kParameterLookupNoWildcard | kParameterLookupNoBaseAtomType;
 768         if (strncmp(com, "bond", len) == 0) {
 769                 BondPar *bp;
 770                 if (sscanf(buf, " %11s %4s %4s %f %f %n", com2, type[0], type[1], &val[0], &val[1], &n) < 5) {
 771                         s_AppendWarning(wbufp, "%s:%d: missing parameter in BOND record\n", fname, lineNumber);
 772                         return 1;
 773                 }
 774                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 775                 itype[1] = AtomTypeEncodeToUInt(type[1]);
 776                 if (itype[0] > itype[1]) {
 777                         i = itype[0];
 778                         itype[0] = itype[1];
 779                         itype[1] = i;
 780                 }
 781                 val[0] *= KCAL2INTERNAL;
 782                 bp = ParameterLookupBondPar(par, itype[0], itype[1], options);
 783                 if (bp != NULL) {
 784                         if (bp->k != val[0] || bp->r0 != val[1]) {
 785                                 s_AppendWarning(wbufp, "%s:%d: The BOND %s-%s parameter appeared twice; the values (%f, %f) are used\n", fname, lineNumber, AtomTypeDecodeToString(itype[0], type[0]), AtomTypeDecodeToString(itype[1], type[1]), val[0], val[1]);
 786                                 retval = 1;
 787                         }
 788                 }
 789                 bp = AssignArray(&par->bondPars, &par->nbondPars, sizeof(*bp), par->nbondPars, NULL);
 790                 bp->type1 = itype[0];
 791                 bp->type2 = itype[1];
 792                 bp->k = val[0];
 793                 bp->r0 = val[1];
 794                 s_StoreComment(kBondParType, bp, buf + n, src);
 795         } else if (strncmp(com, "angle", len) == 0) {
 796                 AnglePar *ap;
 797                 if (sscanf(buf, " %11s %4s %4s %4s %f %f %n", com2, type[0], type[1], type[2], &val[0], &val[1], &n) < 6) {
 798                         s_AppendWarning(wbufp, "%s:%d: missing parameter in ANGLE record\n", fname, lineNumber);
 799                         return 1;
 800                 }
 801                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 802                 itype[1] = AtomTypeEncodeToUInt(type[1]);
 803                 itype[2] = AtomTypeEncodeToUInt(type[2]);
 804                 if (itype[0] > itype[2]) {
 805                         i = itype[0];
 806                         itype[0] = itype[2];
 807                         itype[2] = i;
 808                 }
 809                 val[0] *= KCAL2INTERNAL;
 810                 val[1] *= (3.14159265358979 / 180.0);
 811                 ap = ParameterLookupAnglePar(par, itype[0], itype[1], itype[2], options);
 812                 if (ap != NULL) {
 813                         if (ap->k != val[0] || ap->a0 != val[1]) {
 814                                 s_AppendWarning(wbufp, "%s:%d: The ANGLE %s-%s-%s parameter appeared twice; the values (%f, %f) are used\n", fname, lineNumber, AtomTypeDecodeToString(itype[0], type[0]), AtomTypeDecodeToString(itype[1], type[1]), AtomTypeDecodeToString(itype[2], type[2]), val[0], val[1]);
 815                                 retval = 1;
 816                         }
 817                 }
 818                 ap = AssignArray(&par->anglePars, &par->nanglePars, sizeof(*ap), par->nanglePars, NULL);
 819                 ap->type1 = itype[0];
 820                 ap->type2 = itype[1];
 821                 ap->type3 = itype[2];
 822                 ap->k = val[0];
 823                 ap->a0 = val[1];
 824                 s_StoreComment(kAngleParType, (BondPar *)ap, buf + n, src);
 825         } else if (strncmp(com, "dihedral", len) == 0) {
 826                 TorsionPar *dp;
 827                 if (sscanf(buf, " %11s %4s %4s %4s %4s %f %d %f %n", com2, type[0], type[1], type[2], type[3], &val[0], &ival[0], &val[1], &n) < 8) {
 828                         s_AppendWarning(wbufp, "%s:%d: missing parameter in DIHEDRAL record\n", fname, lineNumber);
 829                         return 1;
 830                 }
 831                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 832                 itype[1] = AtomTypeEncodeToUInt(type[1]);
 833                 itype[2] = AtomTypeEncodeToUInt(type[2]);
 834                 itype[3] = AtomTypeEncodeToUInt(type[3]);
 835                 if (itype[0] > itype[3]) {
 836                         i = itype[0];
 837                         itype[0] = itype[3];
 838                         itype[3] = i;
 839                         i = itype[1];
 840                         itype[1] = itype[2];
 841                         itype[2] = i;
 842                 }
 843                 dp = ParameterLookupDihedralPar(par, itype[0], itype[1], itype[2], itype[3], options);
 844                 val[0] *= KCAL2INTERNAL;
 845                 val[1] *= 3.14159265358979 / 180.0;
 846                 if (dp != NULL) {
 847                         if (dp->mult != 1 || dp->k[0] != val[0] || dp->period[0] != ival[0] || dp->phi0[0] != val[1]) {
 848                                 s_AppendWarning(wbufp, "%s:%d: The DIHEDRAL %s-%s-%s-%s parameter appeared twice; the values (%f, %d, %f) are used\n", fname, lineNumber, AtomTypeDecodeToString(itype[0], type[0]), AtomTypeDecodeToString(itype[1], type[1]), AtomTypeDecodeToString(itype[2], type[2]), AtomTypeDecodeToString(itype[3], type[3]), val[0], ival[0], val[1]);
 849                                 retval = 1;
 850                         }
 851                 }
 852                 dp = AssignArray(&par->dihedralPars, &par->ndihedralPars, sizeof(*dp), par->ndihedralPars, NULL);
 853                 dp->type1 = itype[0];
 854                 dp->type2 = itype[1];
 855                 dp->type3 = itype[2];
 856                 dp->type4 = itype[3];
 857                 dp->k[0] = val[0];
 858                 dp->period[0] = ival[0];
 859                 dp->phi0[0] = val[1];
 860                 dp->mult = 1;
 861                 s_StoreComment(kDihedralParType, (BondPar *)dp, buf + n, src);
 862         } else if (strncmp(com, "improper", len) == 0) {
 863                 TorsionPar *ip;
 864                 if (sscanf(buf, " %11s %4s %4s %4s %4s %f %d %f %n", com2, type[0], type[1], type[2], type[3], &val[0], &ival[0], &val[1], &n) < 8) {
 865                         s_AppendWarning(wbufp, "%s:%d: missing parameter in IMPROPER record\n", fname, lineNumber);
 866                         return 1;
 867                 }
 868                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 869                 itype[1] = AtomTypeEncodeToUInt(type[1]);
 870                 itype[2] = AtomTypeEncodeToUInt(type[2]);
 871                 itype[3] = AtomTypeEncodeToUInt(type[3]);
 872                 if (itype[0] > itype[1]) {
 873                         i = itype[0];
 874                         itype[0] = itype[1];
 875                         itype[1] = i;
 876                 }
 877                 if (itype[0] > itype[3]) {
 878                         i = itype[0];
 879                         itype[0] = itype[3];
 880                         itype[3] = i;
 881                 }
 882                 if (itype[1] > itype[3]) {
 883                         i = itype[1];
 884                         itype[1] = itype[3];
 885                         itype[3] = i;
 886                 }
 887                 ip = ParameterLookupImproperPar(par, itype[0], itype[1], itype[2], itype[3], options);
 888                 val[0] *= KCAL2INTERNAL;
 889                 val[1] *= 3.14159265358979 / 180.0;
 890                 if (ip != NULL) {
 891                         if (ip->mult != 1 || ip->k[0] != val[0] || ip->period[0] != ival[0] || ip->phi0[0] != val[1]) {
 892                                 s_AppendWarning(wbufp, "%s:%d: The IMPROPER %s-%s-%s-%s parameter appeared twice; the values (%f, %d, %f) are used\n", fname, lineNumber, AtomTypeDecodeToString(itype[0], type[0]), AtomTypeDecodeToString(itype[1], type[1]), AtomTypeDecodeToString(itype[2], type[2]), AtomTypeDecodeToString(itype[3], type[3]), val[0], ival[0], val[1]);
 893                                 retval = 1;
 894                         }
 895                 }
 896                 ip = AssignArray(&par->improperPars, &par->nimproperPars, sizeof(*ip), par->nimproperPars, NULL);
 897                 ip->type1 = itype[0];
 898                 ip->type2 = itype[1];
 899                 ip->type3 = itype[2];
 900                 ip->type4 = itype[3];
 901                 ip->k[0] = val[0];
 902                 ip->period[0] = ival[0];
 903                 ip->phi0[0] = val[1];
 904                 ip->mult = 1;
 905                 s_StoreComment(kImproperParType, (BondPar *)ip, buf + n, src);
 906         } else if (strncmp(com, "nonbonded", len) == 0 || strncmp(com, "vdw", len) == 0) {
 907                 VdwPar *vp, vtemp;
 908                 char *p;
 909                 /*  NOTE: the nonbonded record lists "2*sigma", not "sigma"!  */
 910                 int flag = (com[0] == 'v');
 911                 if (sscanf(buf, " %11s %4s %f %f %f %f %n", com2, type[0], &val[0], &val[1], &val[2], &val[3], &n) < 6) {
 912                         s_AppendWarning(wbufp, "%s:%d: missing parameter in %s record\n", fname, lineNumber, (flag ? "VDW" : "NONBONDED"));
 913                         return 1;
 914                 }
 915                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 916                 memset(&vtemp, 0, sizeof(vtemp));
 917                 vtemp.type1 = itype[0];
 918                 vtemp.atomicNumber = 0;  /*  No definition given  */
 919                 vtemp.eps = val[0] * KCAL2INTERNAL;
 920                 vtemp.r_eq = val[1] * (flag ? 1.0 : 0.561231024154687); /* 1/2 * 2**(1/6)  */
 921                 vtemp.A = pow(vtemp.r_eq * 2, 12) * vtemp.eps;
 922                 vtemp.B = 2 * pow(vtemp.r_eq * 2, 6) * vtemp.eps;
 923                 vtemp.eps14 = val[2] * KCAL2INTERNAL;
 924                 vtemp.r_eq14 = val[3] * (flag ? 1.0 : 0.561231024154687); /* 1/2 * 2**(1/6)  */
 925                 vtemp.A14 = pow(vtemp.r_eq14 * 2, 12) * vtemp.eps14;
 926                 vtemp.B14 = 2 * pow(vtemp.r_eq14 * 2, 6) * vtemp.eps14;
 927                 p = buf + n;
 928                 ival[0] = 0;
 929                 val[4] = val[5] = 0.0;
 930                 if (sscanf(p, "%d %n", &ival[0], &n) == 1) {
 931                         vtemp.atomicNumber = ival[0];
 932                         p += n;
 933                 }
 934                 if (sscanf(p, "%f %n", &val[4], &n) == 1) {
 935                         vtemp.weight = val[4];
 936                         p += n;
 937                 }
 938                 if (val[4] == 0.0 && ival[0] != 0)
 939                         vtemp.weight = WeightForAtomicNumber(ival[0]);
 940                 if (sscanf(p, "%f %n", &val[5], &n) == 1) {
 941                         vtemp.polarizability = val[5];
 942                         p += n;
 943                 }
 944                 n = p - buf;
 945                 if (ival[0] == 0 && val[4] != 0.0) {
 946                         for (i = 1; (val[5] = WeightForAtomicNumber(i)) != 0.0; i++) {
 947                                 if (fabs(val[4] - val[5]) < 0.1) {
 948                                         vtemp.atomicNumber = i;
 949                                         break;
 950                                 }
 951                         }
 952                 }
 953                 vp = NULL;
 954                 for (i = 0; i < par->nvdwPars; i++) {
 955                         if (itype[0] == par->vdwPars[i].type1) {
 956                                 vp = par->vdwPars + i;
 957                                 if (vp->A != vtemp.A || vp->B != vtemp.B || vp->A14 != vtemp.A14 || vp->B14 != vtemp.B14) {
 958                                         s_AppendWarning(wbufp, "%s:%d: The %s %s parameter appeared twice; the values (%f, %f, %f, %f, %d, %f, %f) are used\n", fname, lineNumber, (flag ? "VDW" : "NONBONDED"), AtomTypeDecodeToString(itype[0], type[0]), val[0], val[1], val[2], val[3], ival[0], val[4], val[5]);
 959                                         retval = 1;
 960                                 }
 961                                 break;
 962                         }
 963                 }
 964                 vp = AssignArray(&par->vdwPars, &par->nvdwPars, sizeof(*vp), i, NULL);
 965                 vtemp.com = vp->com;  /*  Keep comment field  */
 966                 *vp = vtemp;
 967                 s_StoreComment(kVdwParType, (BondPar *)vp, buf + n, src);
 968         } else if (strncmp(com, "nbfi", len) == 0 || strncmp(com, "vdwpair", len) == 0) {
 969                 VdwPairPar *vp, vtemp;
 970                 int flag = (com[0] == 'v');
 971                 if (sscanf(buf, " %11s %4s %4s %f %f %f %f %n", com2, type[0], type[1], &val[0], &val[1], &val[2], &val[3], &n) < 6) {
 972                         s_AppendWarning(wbufp, "%s:%d: missing parameter in %s record\n", fname, lineNumber, (flag ? "VDWP" : "NBFI"));
 973                         return 1;
 974                 }
 975                 itype[0] = AtomTypeEncodeToUInt(type[0]);
 976                 itype[1] = AtomTypeEncodeToUInt(type[1]);
 977                 if (itype[0] > itype[1]) {
 978                         i = itype[0];
 979                         itype[0] = itype[1];
 980                         itype[1] = i;
 981                 }
 982                 vtemp.type1 = itype[0];
 983                 vtemp.type2 = itype[1];
 984                 if (flag) {  /*  eps/r_eq representation  */
 985                         vtemp.eps = val[0] * KCAL2INTERNAL;
 986                         vtemp.r_eq = val[1];
 987                         vtemp.A = pow(val[1] * 2, 12) * vtemp.eps;
 988                         vtemp.B = 2 * pow(val[1] * 2, 6) * vtemp.eps;
 989                         vtemp.eps14 = val[2] * KCAL2INTERNAL;
 990                         vtemp.r_eq14 = val[3];
 991                         vtemp.A14 = pow(val[3] * 2, 12) * vtemp.eps14;
 992                         vtemp.B14 = 2 * pow(val[3] * 2, 6) * vtemp.eps14;
 993                 } else {    /*  A/B representation  */
 994                         vtemp.A = val[0] * KCAL2INTERNAL;
 995                         vtemp.B = val[1] * KCAL2INTERNAL;
 996                         vtemp.eps = pow(0.25 * vtemp.B * vtemp.B / vtemp.A, 0.16666666667);
 997                         vtemp.r_eq = pow(vtemp.A / vtemp.B * 2.0, 0.16666666667) * 0.5;
 998                         vtemp.A14 = val[2] * KCAL2INTERNAL;
 999                         vtemp.B14 = val[3] * KCAL2INTERNAL;
1000                         vtemp.eps14 = pow(0.25 * vtemp.B14 * vtemp.B14 / vtemp.A14, 0.16666666667);
1001                         vtemp.r_eq14 = pow(vtemp.A14 / vtemp.B14 * 2.0, 0.16666666667) * 0.5;
1002                 }
1003                 vp = NULL;
1004                 for (i = 0; i < par->nvdwpPars; i++) {
1005                         if (itype[0] == par->vdwpPars[i].type1 && itype[1] == par->vdwpPars[i].type2) {
1006                                 vp = par->vdwpPars + i;
1007                                 if (vp->A != vtemp.A || vp->B != vtemp.B || vp->A14 != vtemp.A14 || vp->B14 != vtemp.B14) {
1008                                         s_AppendWarning(wbufp, "%s:%d: The %s %s-%s parameter appeared twice; the values (%f, %f, %f, %f) are used\n", fname, lineNumber, (flag ? "VDWP" : "NBFI"), AtomTypeDecodeToString(itype[0], type[0]), AtomTypeDecodeToString(itype[1], type[1]), val[0], val[1], val[2], val[3]);
1009                                         retval = 1;
1010                                 }
1011                                 break;
1012                         }
1013                 }
1014                 vp = AssignArray(&par->vdwpPars, &par->nvdwpPars, sizeof(*vp), i, NULL);
1015                 vtemp.com = vp->com;  /*  Keep comment field  */
1016                 *vp = vtemp;
1017                 s_StoreComment(kVdwPairParType, (BondPar *)vp, buf + n, src);
1018         } else {
1019                 s_AppendWarning(wbufp, "%s:%d: unknown keyword %s\n", fname, lineNumber, com);
1020                 return 1;
1021         }
1022         return retval;
1023 }
1024
1025 int
1026 ParameterReadFromFile(Parameter *par, const char *fname, char **outWarningMessage, int *outWarningCount)
1027 {
1028         char *wbuf;
1029         char buf[1024];
1030         FILE *fp = NULL;
1031         int first = 0;
1032         int wcount;
1033         int lineNumber;
1034         int src_idx;
1035
1036         if (par == NULL) {
1037                 par = gBuiltinParameters;
1038                 if (par == NULL) {
1039                         par = ParameterNew();
1040                         gBuiltinParameters = par;
1041                         first = 1;
1042                 }
1043         }
1044
1045         wcount = 0;
1046         wbuf = NULL;
1047
1048         fp = fopen(fname, "rb");
1049         if (fp == NULL) {
1050                 s_AppendWarning(&wbuf, "Cannot open parameter file %s\n", fname);
1051                 wcount = -1;
1052                 goto exit;
1053         }
1054
1055         if (par != gBuiltinParameters || first)
1056                 src_idx = 0;
1057         else
1058                 src_idx = ParameterCommentIndexForGlobalFileName(fname);
1059
1060         if (par == gBuiltinParameters && !first) {
1061                 /*  Ensure the "source" field is unique  */
1062                 int i;
1063                 ParFileInfo *ip;
1064                 const char *p;
1065                 for (i = 0; i < gGlobalParInfo.count; i++) {
1066                         if (gGlobalParInfo.files[i].src == src_idx)
1067                                 break;
1068                 }
1069                 if (i < gGlobalParInfo.count) {
1070                         /*  Delete the existing Parameters from the same source  */
1071                         ParameterDeleteAllEntriesForSource(par, src_idx);
1072                 }
1073
1074                 /*  Register the global file info  */
1075                 ip = AssignArray(&(gGlobalParInfo.files), &(gGlobalParInfo.count), sizeof(ParFileInfo), gGlobalParInfo.count, NULL);
1076                 for (p = fname + strlen(fname) - 1; p >= fname; p--) {
1077                         if (*p == '/' || *p == '\\')
1078                                 break;
1079                 }
1080                 if (p < fname)
1081                         ip->dir = NULL;
1082                 else {
1083                         i = p - fname;
1084                         ip->dir = (char *)malloc(i + 1);
1085                         if (ip->dir != NULL) {
1086                                 strncpy(ip->dir, fname, i);
1087                                 ip->dir[i] = 0;
1088                         }
1089                 }
1090                 p++;
1091                 ip->name = strdup(p);
1092                 ip->src = src_idx;
1093         }
1094
1095         lineNumber = 0;
1096         while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
1097                 int wc1 = ParameterReadFromString(par, buf, &wbuf, fname, lineNumber, src_idx);
1098                 if (wc1 >= 0)
1099                         wcount += wc1;
1100                 else {
1101                         wcount = wc1;
1102                         break;
1103                 }
1104         }
1105         if (first)
1106                 gGlobalParInfo.builtinCount = gGlobalParInfo.count;
1107
1108 exit:
1109         if (fp != NULL)
1110                 fclose(fp);
1111         if (outWarningMessage != NULL)
1112                 *outWarningMessage = wbuf;
1113         else if (wbuf != NULL)
1114                 free(wbuf);
1115         if (outWarningCount != NULL)
1116                 *outWarningCount = (wcount >= 0 ? wcount : 0);
1117
1118         return (wcount >= 0 ? 0 : wcount);
1119 }
1120
1121 int
1122 ParameterAppendToFile(Parameter *par, FILE *fp)
1123 {
1124         int i, n;
1125         char cbuf[6][8];
1126         char buf[120];
1127         int bufsize = sizeof(buf);
1128
1129         if (par == NULL)
1130                 return 0;
1131
1132         n = 0;
1133         if (par->nbondPars > 0)
1134                 fprintf(fp, "! type1 type2 k r0\n");
1135         for (i = 0; i < par->nbondPars; i++) {
1136                 BondPar *bp = par->bondPars + i;
1137                 AtomTypeDecodeToString(bp->type1, cbuf[0]);
1138                 AtomTypeDecodeToString(bp->type2, cbuf[1]);
1139                 s_CommentToString(buf, bufsize, bp);
1140                 fprintf(fp, "bond %s %s %.6f %f%s\n", cbuf[0], cbuf[1], bp->k * INTERNAL2KCAL, bp->r0, buf);
1141                 n++;
1142         }
1143         if (par->nanglePars > 0)
1144                 fprintf(fp, "! type1 type2 type3 k a0\n");
1145         for (i = 0; i < par->nanglePars; i++) {
1146                 AnglePar *ap = par->anglePars + i;
1147                 AtomTypeDecodeToString(ap->type1, cbuf[0]);
1148                 AtomTypeDecodeToString(ap->type2, cbuf[1]);
1149                 AtomTypeDecodeToString(ap->type3, cbuf[2]);
1150                 s_CommentToString(buf, bufsize, ap);
1151                 fprintf(fp, "angle %s %s %s %.6f %f%s\n", cbuf[0], cbuf[1], cbuf[2], ap->k * INTERNAL2KCAL, ap->a0 * kRad2Deg, buf);
1152                 n++;
1153         }
1154         if (par->ndihedralPars > 0)
1155                 fprintf(fp, "! type1 type2 type3 type4 k periodicity phi0\n");
1156         for (i = 0; i < par->ndihedralPars; i++) {
1157                 TorsionPar *tp = par->dihedralPars + i;
1158                 AtomTypeDecodeToString(tp->type1, cbuf[0]);
1159                 AtomTypeDecodeToString(tp->type2, cbuf[1]);
1160                 AtomTypeDecodeToString(tp->type3, cbuf[2]);
1161                 AtomTypeDecodeToString(tp->type4, cbuf[3]);
1162                 s_CommentToString(buf, bufsize, tp);
1163                 fprintf(fp, "dihe %s %s %s %s %.6f %d %f%s\n", cbuf[0], cbuf[1], cbuf[2], cbuf[3], tp->k[0] * INTERNAL2KCAL, tp->period[0], tp->phi0[0] * kRad2Deg, buf);
1164                 n++;
1165         }
1166         if (par->nimproperPars > 0)
1167                 fprintf(fp, "! type1 type2 type3 type4 k periodicity phi0\n");
1168         for (i = 0; i < par->nimproperPars; i++) {
1169                 TorsionPar *tp = par->improperPars + i;
1170                 AtomTypeDecodeToString(tp->type1, cbuf[0]);
1171                 AtomTypeDecodeToString(tp->type2, cbuf[1]);
1172                 AtomTypeDecodeToString(tp->type3, cbuf[2]);
1173                 AtomTypeDecodeToString(tp->type4, cbuf[3]);
1174                 s_CommentToString(buf, bufsize, tp);
1175                 fprintf(fp, "impr %s %s %s %s %.6f %d %f%s\n", cbuf[0], cbuf[1], cbuf[2], cbuf[3], tp->k[0] * INTERNAL2KCAL, tp->period[0], tp->phi0[0] * kRad2Deg, buf);
1176                 n++;
1177         }
1178         if (par->nvdwPars > 0)
1179                 fprintf(fp, "! type eps r_eq eps14 r_eq14 atomic_number weight\n");
1180         for (i = 0; i < par->nvdwPars; i++) {
1181                 VdwPar *vp = par->vdwPars + i;
1182         /*      Double eps, eps14;  */
1183                 AtomTypeDecodeToString(vp->type1, cbuf[0]);
1184         /*      eps = (vp->A == 0.0 ? 0.0 : vp->B * vp->B / vp->A * 0.25 * INTERNAL2KCAL);
1185                 eps14 = (vp->A14 == 0.0 ? 0.0 : vp->B14 * vp->B14 / vp->A14 * 0.25 * INTERNAL2KCAL);  */
1186                 s_CommentToString(buf, bufsize, vp);
1187                 fprintf(fp, "vdw %s %.6f %.6f %.6f %.6f %d %f%s\n", cbuf[0], vp->eps * INTERNAL2KCAL, vp->r_eq, vp->eps14 * INTERNAL2KCAL, vp->r_eq14, vp->atomicNumber, vp->weight, buf);  /*  polarizability is not written because it is not used now  */
1188                 n++;
1189         }
1190         if (par->nvdwpPars > 0)
1191                 fprintf(fp, "! type1 type2 eps r_eq eps14 r_eq14\n");
1192         for (i = 0; i < par->nvdwpPars; i++) {
1193                 VdwPairPar *vpp = par->vdwpPars + i;
1194         /*      Double eps, eps14;  */
1195                 AtomTypeDecodeToString(vpp->type1, cbuf[0]);
1196                 AtomTypeDecodeToString(vpp->type2, cbuf[1]);
1197         /*      eps = (vpp->A == 0.0 ? 0.0 : vpp->B * vpp->B / vpp->A * 0.25 * INTERNAL2KCAL);
1198                 eps14 = (vpp->A14 == 0.0 ? 0.0 : vpp->B14 * vpp->B14 / vpp->A14 * 0.25 * INTERNAL2KCAL); */
1199                 s_CommentToString(buf, bufsize, vpp);
1200                 fprintf(fp, "vdwp %s %s %.6f %.6f %.6f %.6f%s\n", cbuf[0], cbuf[1], vpp->eps * INTERNAL2KCAL, vpp->r_eq, vpp->eps14 * INTERNAL2KCAL, vpp->r_eq14, buf);
1201                 n++;
1202         }
1203 /*      if (par->natomPars > 0)
1204                 fprintf(fp, "! name atomic_number radius red green blue weight\n");
1205         for (i = 0; i < par->natomPars; i++) {
1206                 ElementPar *app = par->atomPars + i;
1207                 s_CommentToString(buf, bufsize, app);
1208                 fprintf(fp, "element %.4s %d %f %f %f %f %f%s\n", app->name, app->number, app->radius, app->r, app->g, app->b, app->weight, buf);
1209                 n++;
1210         } */
1211         return n;
1212 }
1213
1214 int
1215 ParameterDeleteAllEntriesForSource(Parameter *par, int src_idx)
1216 {
1217         int i, n, type;
1218         UnionPar *up;
1219         IntGroup *ig;
1220 /*      if (fname == NULL)
1221                 return -1;
1222         src_idx = ParameterCommentIndexForGlobalFileName(fname); */
1223         if (src_idx == 0)
1224                 return -1;
1225         n = 0;
1226         for (type = kFirstParType; type <= kLastParType; type++) {
1227                 ig = IntGroupNew();
1228                 for (i = ParameterGetCountForType(par, type) - 1; i >= 0; i--) {
1229                         up = ParameterGetUnionParFromTypeAndIndex(par, type, i);
1230                         if (up != NULL && up->bond.src == src_idx)
1231                                 IntGroupAdd(ig, i, 1);
1232                 }
1233                 i = IntGroupGetCount(ig);
1234                 if (i > 0) {
1235                         ParameterDelete(par, type, NULL, ig);
1236                         n += i;
1237                 }
1238                 IntGroupRelease(ig);
1239         }
1240         if (par == gBuiltinParameters) {
1241                 /*  Unregister from the global info  */
1242                 for (i = gGlobalParInfo.builtinCount; i < gGlobalParInfo.count; i++) {
1243                         if (gGlobalParInfo.files[i].src == src_idx) {
1244                                 DeleteArray(&gGlobalParInfo.files, &gGlobalParInfo.count, sizeof(ParFileInfo), i, 1, NULL);
1245                                 break;
1246                         }
1247                 }
1248         }
1249         return n;
1250 }
1251
1252 #pragma mark ====== Parameter Comments ======
1253
1254 static const char **sParameterComments;
1255 static Int sNumberOfParameterComments;
1256
1257 int
1258 ParameterCommentIndex(const char *comment)
1259 {
1260         int i;
1261         char *p, *pp;
1262         if (comment == NULL || comment[0] == 0)
1263                 return 0;
1264         /*  Duplicate the comment, convert whitespaces to ' ', and chop trailing whitespaces  */
1265         p = strdup(comment);
1266         i = 0;
1267         for (pp = p + strlen(p) - 1; pp >= p; pp--) {
1268                 if (isspace(*pp)) {
1269                         if (i == 0)
1270                                 *pp = 0;
1271                         else *pp = ' ';
1272                 } else i++;
1273         }
1274         for (i = 1; i < sNumberOfParameterComments; i++) {
1275                 if (strcmp(sParameterComments[i], p) == 0) {
1276                         free(p);
1277                         return i;
1278                 }
1279         }
1280         if (sNumberOfParameterComments == 0) {
1281                 /*  Index 0 is skipped  */
1282                 AssignArray(&sParameterComments, &sNumberOfParameterComments, sizeof(char *), 1, &p);
1283         } else {
1284                 AssignArray(&sParameterComments, &sNumberOfParameterComments, sizeof(char *), sNumberOfParameterComments, &p);
1285         }
1286         return sNumberOfParameterComments - 1;
1287 }
1288
1289 const char *
1290 ParameterGetComment(int idx)
1291 {
1292         if (idx <= 0 || idx >= sNumberOfParameterComments)
1293                 return NULL;  /*  No such number  */
1294         return sParameterComments[idx];
1295 }
1296
1297 #pragma mark ====== Parameter Lookup ======
1298
1299 #define s_ParMatch(_t1, _t2, _nowildcard) (_t1 == _t2 || (!_nowildcard && _t1 == kAtomTypeWildcard))
1300
1301 /*  Returns non-zero if the parameter record contains designated atom_type.
1302  The atom_type can also be an atom index.  */
1303 int
1304 ParameterDoesContainAtom(Int type, UnionPar *up, UInt atom_type, Int options)
1305 {
1306 #define CHECK_FIELD(_tp, _field) s_ParMatch((((_tp *)up)->_field), atom_type, nowildcard)
1307         Int nowildcard = (options & kParameterLookupNoWildcard);
1308         switch (type) {
1309                 case kBondParType:
1310                         return CHECK_FIELD(BondPar, type1) || CHECK_FIELD(BondPar, type2);
1311                 case kAngleParType:
1312                         return CHECK_FIELD(AnglePar, type1) || CHECK_FIELD(AnglePar, type2) || CHECK_FIELD(AnglePar, type3);
1313                 case kDihedralParType:
1314                 case kImproperParType:
1315                         return CHECK_FIELD(TorsionPar, type1) || CHECK_FIELD(TorsionPar, type2) || CHECK_FIELD(TorsionPar, type3) || CHECK_FIELD(TorsionPar, type4);
1316                 case kVdwParType:
1317                         return CHECK_FIELD(VdwPar, type1);
1318                 case kVdwPairParType:
1319                         return CHECK_FIELD(VdwPairPar, type1) || CHECK_FIELD(VdwPairPar, type2);
1320                 case kVdwCutoffParType:
1321                         if (((VdwCutoffPar *)up)->type == 1)
1322                                 return CHECK_FIELD(VdwCutoffPar, n1) || CHECK_FIELD(VdwCutoffPar, n2) || CHECK_FIELD(VdwCutoffPar, n3) || CHECK_FIELD(VdwCutoffPar, n4);
1323                         else return 0;
1324                 default: return 0;
1325         }
1326 }
1327
1328 BondPar *
1329 ParameterLookupBondPar(Parameter *par, UInt t1, UInt t2, Int options)
1330 {
1331         int i;
1332         BondPar *bp;
1333         Int nowildcard = (options & kParameterLookupNoWildcard);
1334         if (par == NULL)
1335                 return NULL;
1336         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1337                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1338         for (i = par->nbondPars - 1, bp = par->bondPars + i; i >= 0; i--, bp--) {
1339                 if ((bp->src > 0 && !(options & kParameterLookupGlobal))
1340                         || (bp->src == 0 && !(options & kParameterLookupLocal))
1341                         || (bp->src < 0 && !(options & kParameterLookupMissing)))
1342                         continue;
1343                 if (s_ParMatch(bp->type1, t1, nowildcard) && s_ParMatch(bp->type2, t2, nowildcard))
1344                         return bp;
1345                 if (s_ParMatch(bp->type1, t2, nowildcard) && s_ParMatch(bp->type2, t1, nowildcard))
1346                         return bp;
1347         }
1348         if (options & kParameterLookupNoBaseAtomType)
1349                 return NULL;
1350         return ParameterLookupBondPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1351 }
1352
1353 AnglePar *
1354 ParameterLookupAnglePar(Parameter *par, UInt t1, UInt t2, UInt t3, Int options)
1355 {
1356         int i;
1357         AnglePar *ap;
1358         Int nowildcard = (options & kParameterLookupNoWildcard);
1359         if (par == NULL)
1360                 return NULL;
1361         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1362                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1363         for (i = par->nanglePars - 1, ap = par->anglePars + i; i >= 0; i--, ap--) {
1364                 if ((ap->src > 0 && !(options & kParameterLookupGlobal))
1365                         || (ap->src == 0 && !(options & kParameterLookupLocal))
1366                         || (ap->src < 0 && !(options & kParameterLookupMissing)))
1367                         continue;
1368                 if (s_ParMatch(ap->type1, t1, nowildcard) && s_ParMatch(ap->type2, t2, nowildcard) && s_ParMatch(ap->type3, t3, nowildcard))
1369                         return ap;
1370                 if (s_ParMatch(ap->type1, t3, nowildcard) && s_ParMatch(ap->type2, t2, nowildcard) && s_ParMatch(ap->type3, t1, nowildcard))
1371                         return ap;
1372         }
1373         if (options & kParameterLookupNoBaseAtomType)
1374                 return NULL;
1375         return ParameterLookupAnglePar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1376 }
1377
1378 TorsionPar *
1379 ParameterLookupDihedralPar(Parameter *par, UInt t1, UInt t2, UInt t3, UInt t4, Int options)
1380 {
1381         int i;
1382         TorsionPar *tp;
1383         Int nowildcard = (options & kParameterLookupNoWildcard);
1384         if (par == NULL)
1385                 return NULL;
1386         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1387                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1388         for (i = par->ndihedralPars - 1, tp = par->dihedralPars + i; i >= 0; i--, tp--) {
1389                 if ((tp->src > 0 && !(options & kParameterLookupGlobal))
1390                         || (tp->src == 0 && !(options & kParameterLookupLocal))
1391                         || (tp->src < 0 && !(options & kParameterLookupMissing)))
1392                         continue;
1393                 if (s_ParMatch(tp->type1, t1, nowildcard) && s_ParMatch(tp->type2, t2, nowildcard) && s_ParMatch(tp->type3, t3, nowildcard) && s_ParMatch(tp->type4, t4, nowildcard))
1394                         return tp;
1395                 if (s_ParMatch(tp->type1, t4, nowildcard) && s_ParMatch(tp->type2, t3, nowildcard) && s_ParMatch(tp->type3, t2, nowildcard) && s_ParMatch(tp->type4, t1, nowildcard))
1396                         return tp;
1397         }
1398         if (options & kParameterLookupNoBaseAtomType)
1399                 return NULL;
1400         return ParameterLookupDihedralPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, t4 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1401 }
1402
1403 TorsionPar *
1404 ParameterLookupImproperPar(Parameter *par, UInt t1, UInt t2, UInt t3, UInt t4, Int options)
1405 {
1406         int i;
1407         TorsionPar *tp;
1408         Int nowildcard = (options & kParameterLookupNoWildcard);
1409         if (par == NULL)
1410                 return NULL;
1411         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1412                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1413         for (i = par->nimproperPars - 1, tp = par->improperPars + i; i >= 0; i--, tp--) {
1414                 if ((tp->src > 0 && !(options & kParameterLookupGlobal))
1415                         || (tp->src == 0 && !(options & kParameterLookupLocal))
1416                         || (tp->src < 0 && !(options & kParameterLookupMissing)))
1417                         continue;
1418                 if (!s_ParMatch(tp->type3, t3, nowildcard))
1419                         continue;
1420                 if ((s_ParMatch(tp->type1, t1, nowildcard) && s_ParMatch(tp->type2, t2, nowildcard) && s_ParMatch(tp->type4, t4, nowildcard))
1421                         || (s_ParMatch(tp->type1, t1, nowildcard) && s_ParMatch(tp->type2, t4, nowildcard) && s_ParMatch(tp->type4, t2, nowildcard))
1422                         || (s_ParMatch(tp->type1, t2, nowildcard) && s_ParMatch(tp->type2, t1, nowildcard) && s_ParMatch(tp->type4, t4, nowildcard))
1423                         || (s_ParMatch(tp->type1, t2, nowildcard) && s_ParMatch(tp->type2, t4, nowildcard) && s_ParMatch(tp->type4, t1, nowildcard))
1424                         || (s_ParMatch(tp->type1, t4, nowildcard) && s_ParMatch(tp->type2, t1, nowildcard) && s_ParMatch(tp->type4, t2, nowildcard))
1425                         || (s_ParMatch(tp->type1, t4, nowildcard) && s_ParMatch(tp->type2, t2, nowildcard) && s_ParMatch(tp->type4, t1, nowildcard)))
1426                         return tp;
1427         }
1428         if (options & kParameterLookupNoBaseAtomType)
1429                 return NULL;
1430         return ParameterLookupImproperPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, t4 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1431 }
1432
1433 VdwPar *
1434 ParameterLookupVdwPar(Parameter *par, UInt t1, Int options)
1435 {
1436         int i;
1437         VdwPar *vp;
1438         Int nowildcard = (options & kParameterLookupNoWildcard);
1439         if (par == NULL)
1440                 return NULL;
1441         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1442                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1443         for (i = par->nvdwPars - 1, vp = par->vdwPars + i; i >= 0; i--, vp--) {
1444                 if ((vp->src > 0 && !(options & kParameterLookupGlobal))
1445                         || (vp->src == 0 && !(options & kParameterLookupLocal))
1446                         || (vp->src < 0 && !(options & kParameterLookupMissing)))
1447                         continue;
1448                 if (s_ParMatch(vp->type1, t1, nowildcard))
1449                         return vp;
1450         }
1451         if (options & kParameterLookupNoBaseAtomType)
1452                 return NULL;
1453         return ParameterLookupVdwPar(par, t1 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1454 }
1455
1456 VdwPairPar *
1457 ParameterLookupVdwPairPar(Parameter *par, UInt t1, UInt t2, Int options)
1458 {
1459         int i;
1460         VdwPairPar *vp;
1461         Int nowildcard = (options & kParameterLookupNoWildcard);
1462         if (par == NULL)
1463                 return NULL;
1464         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1465                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1466         for (i = par->nvdwpPars - 1, vp = par->vdwpPars + i; i >= 0; i--, vp--) {
1467                 if ((vp->src > 0 && !(options & kParameterLookupGlobal))
1468                         || (vp->src == 0 && !(options & kParameterLookupLocal))
1469                         || (vp->src < 0 && !(options & kParameterLookupMissing)))
1470                         continue;
1471                 if ((s_ParMatch(vp->type1, t1, nowildcard) && s_ParMatch(vp->type2, t2, nowildcard))
1472                         || (s_ParMatch(vp->type1, t2, nowildcard) && s_ParMatch(vp->type2, t1, nowildcard)))
1473                         return vp;
1474         }
1475         if (options & kParameterLookupNoBaseAtomType)
1476                 return NULL;
1477         return ParameterLookupVdwPairPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1478 }
1479
1480 VdwCutoffPar *
1481 ParameterLookupVdwCutoffPar(Parameter *par, UInt t1, UInt t2, Int options)
1482 {
1483         int i;
1484         VdwCutoffPar *vp;
1485         Int nowildcard = (options & kParameterLookupNoWildcard);
1486         if (par == NULL)
1487                 return NULL;
1488         if ((options & (kParameterLookupGlobal | kParameterLookupLocal | kParameterLookupMissing)) == 0)
1489                 options |= kParameterLookupGlobal | kParameterLookupLocal;
1490         for (i = par->nvdwCutoffPars - 1, vp = par->vdwCutoffPars + i; i >= 0; i--, vp--) {
1491                 if ((vp->src > 0 && !(options & kParameterLookupGlobal))
1492                         || (vp->src == 0 && !(options & kParameterLookupLocal))
1493                         || (vp->src < 0 && !(options & kParameterLookupMissing)))
1494                         continue;
1495                 if (vp->type == 0) {
1496                         if (s_ParMatch(vp->n1, t1, nowildcard) && s_ParMatch(vp->n2, t2, nowildcard))
1497                                 return vp;
1498                         if (s_ParMatch(vp->n1, t2, nowildcard) && s_ParMatch(vp->n2, t1, nowildcard))
1499                                 return vp;
1500                 } else {
1501                         if (vp->n1 <= t1 && vp->n2 >= t1 && vp->n3 <= t2 && vp->n4 <= t2)
1502                                 return vp;
1503                         if (vp->n1 <= t2 && vp->n2 >= t2 && vp->n3 <= t1 && vp->n4 >= t1)
1504                                 return vp;
1505                 }
1506         }
1507         if (options & kParameterLookupNoBaseAtomType)
1508                 return NULL;
1509         return ParameterLookupVdwCutoffPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1510 }
1511
1512 #pragma mark ====== Table View Support ======
1513
1514 int
1515 ParameterTableNumberOfRows(Parameter *par)
1516 {
1517         if (par == NULL)
1518                 return 0;
1519         return par->nvdwPars + par->nbondPars + par->nanglePars + par->ndihedralPars + par->nimproperPars + par->nvdwpPars + 6;
1520 }
1521
1522 int
1523 ParameterTableGetItemIndex(Parameter *par, int row, int *type)
1524 {
1525         if (par == NULL || row < 0) {
1526                 *type = kInvalidParType;
1527                 return -1;
1528         }
1529         if (--row < par->nvdwPars) {
1530                 *type = kVdwParType;
1531         } else if ((row -= par->nvdwPars + 1) < par->nbondPars) {
1532                 *type = kBondParType;
1533         } else if ((row -= par->nbondPars + 1) < par->nanglePars) {
1534                 *type = kAngleParType;
1535         } else if ((row -= par->nanglePars + 1) < par->ndihedralPars) {
1536                 *type = kDihedralParType;
1537         } else if ((row -= par->ndihedralPars + 1) < par->nimproperPars) {
1538                 *type = kImproperParType;
1539         } else if ((row -= par->nimproperPars + 1) < par->nvdwpPars) {
1540                 *type = kVdwPairParType;
1541         } else {
1542                 *type = kInvalidParType;
1543                 return -1;
1544         }
1545         return row;
1546 }
1547
1548 int
1549 ParameterTableGetRowFromTypeAndIndex(Parameter *par, int type, int idx)
1550 {
1551         int ofs = 1;
1552         if (idx < -1)
1553                 return -1;  /*  Invalid  */
1554         if (type == kVdwParType)
1555                 return (idx < par->nvdwPars ? idx + ofs : -1);
1556         ofs += par->nvdwPars + 1;
1557         if (type == kBondParType)
1558                 return (idx < par->nbondPars ? idx + ofs : -1);
1559         ofs += par->nbondPars + 1;
1560         if (type == kAngleParType)
1561                 return (idx < par->nanglePars ? idx + ofs : -1);
1562         ofs += par->nanglePars + 1;
1563         if (type == kDihedralParType)
1564                 return (idx < par->ndihedralPars ? idx + ofs : -1);
1565         ofs += par->ndihedralPars + 1;
1566         if (type == kImproperParType)
1567                 return (idx < par->nimproperPars ? idx + ofs : -1);
1568         ofs += par->nimproperPars + 1;
1569         if (type == kVdwPairParType)
1570                 return (idx < par->nvdwpPars ? idx + ofs : -1);
1571         return -1;
1572 }
1573
1574 UnionPar *
1575 ParameterTableGetItemPtr(Parameter *par, int row, int *type)
1576 {
1577         if (par == NULL || row < 0) {
1578                 *type = kInvalidParType;
1579                 return NULL;
1580         }
1581         if (--row < par->nvdwPars) {
1582                 *type = kVdwParType;
1583                 return (UnionPar *)(row >= 0 ? par->vdwPars + row : NULL);
1584         } else if ((row -= par->nvdwPars + 1) < par->nbondPars) {
1585                 *type = kBondParType;
1586                 return (UnionPar *)(row >= 0 ? par->bondPars + row : NULL);
1587         } else if ((row -= par->nbondPars + 1) < par->nanglePars) {
1588                 *type = kAngleParType;
1589                 return (UnionPar *)(row >= 0 ? par->anglePars + row : NULL);
1590         } else if ((row -= par->nanglePars + 1) < par->ndihedralPars) {
1591                 *type = kDihedralParType;
1592                 return (UnionPar *)(row >= 0 ? par->dihedralPars + row : NULL);
1593         } else if ((row -= par->ndihedralPars + 1) < par->nimproperPars) {
1594                 *type = kImproperParType;
1595                 return (UnionPar *)(row >= 0 ? par->improperPars + row : NULL);
1596         } else if ((row -= par->nimproperPars + 1) < par->nvdwpPars) {
1597                 *type = kVdwPairParType;
1598                 return (UnionPar *)(row >= 0 ? par->vdwpPars + row : NULL);
1599 /*      } else if ((row -= par->nvdwpPars + 1) < par->natomPars) {
1600                 *type = kElementParType;
1601                 return (UnionPar *)(row >= 0 ? par->atomPars + row : NULL); */
1602         } else {
1603                 *type = kInvalidParType;
1604                 return NULL;
1605         }
1606 }
1607
1608 void
1609 ParameterTableGetItemText(Parameter *par, int column, int row, char *buf, int bufsize)
1610 {
1611         static char *sBondParTitles[] = {"", "Bonds", "k", "r0"};
1612         static char *sAngleParTitles[] = {"", "Angles", "k", "a0"};
1613         static char *sDihedralParTitles[] = {"", "Dihedrals", "k", "period", "phi0"};
1614         static char *sImproperParTitles[] = {"", "Impropers", "k", "period", "phi0"};
1615         static char *sVdwParTitles[] = {"", "VDWs", "eps", "r", "eps14", "r14", "atomNo", "weight"};
1616         static char *sVdwPairParTitles[] = {"", "VDW Pairs", "eps", "r", "eps14", "r14"};
1617 /*      static char *sAtomParTitles[] = {"", "Atom Display", "atomNo", "radius", "red", "green", "blue", "weight"}; */
1618         const char *p;
1619         int type;
1620         UnionPar *up = NULL;
1621         char types[4][8];
1622         buf[0] = 0;
1623         if (par == NULL || row < 0)
1624                 return;
1625         up = ParameterTableGetItemPtr(par, row, &type);
1626         switch (type) {
1627                 case kVdwParType: {
1628                         VdwPar *vp = (VdwPar *)up;
1629                         if (vp == NULL) {
1630                                 if (column >= 0 && column < 8)
1631                                         snprintf(buf, bufsize, "%s", sVdwParTitles[column]);
1632                                 return;
1633                         }
1634                         switch (column) {
1635                                 case 0: snprintf(buf, bufsize, "vdw"); break;
1636                                 case 1:
1637                                         AtomTypeDecodeToString(vp->type1, types[0]);
1638                                         snprintf(buf, bufsize, "%s", types[0]);
1639                                         break;
1640                                 case 2:
1641                                         snprintf(buf, bufsize, "%.5f", vp->eps * INTERNAL2KCAL);
1642                                         break;
1643                                 case 3:
1644                                         snprintf(buf, bufsize, "%.5f", vp->r_eq);
1645                                         break;
1646                                 case 4:
1647                                         snprintf(buf, bufsize, "%.5f", vp->eps14 * INTERNAL2KCAL);
1648                                         break;
1649                                 case 5:
1650                                         snprintf(buf, bufsize, "%.5f", vp->r_eq14);
1651                                         break;
1652                                 case 6:
1653                                         snprintf(buf, bufsize, "%d", vp->atomicNumber);
1654                                         break;
1655                                 case 7:
1656                                         snprintf(buf, bufsize, "%.3f", vp->weight);
1657                                         break;
1658                         }
1659                         break;
1660                 }
1661                 case kBondParType: {
1662                         BondPar *bp = (BondPar *)up;
1663                         if (bp == NULL) {
1664                                 if (column >= 0 && column < 4)
1665                                         snprintf(buf, bufsize, "%s", sBondParTitles[column]);
1666                                 return;
1667                         }
1668                         switch (column) {
1669                                 case 0: snprintf(buf, bufsize, "bond"); break;
1670                                 case 1:
1671                                         AtomTypeDecodeToString(bp->type1, types[0]);
1672                                         AtomTypeDecodeToString(bp->type2, types[1]);
1673                                         snprintf(buf, bufsize, "%s-%s", types[0], types[1]);
1674                                         break;
1675                                 case 2:
1676                                 case 3:
1677                                         snprintf(buf, bufsize, "%.3f", (column == 2 ? bp->k * INTERNAL2KCAL : bp->r0));
1678                                         break;
1679                         }
1680                         break;
1681                 }
1682                 case kAngleParType: {
1683                         AnglePar *ap = (AnglePar *)up;
1684                         if (ap == NULL) {
1685                                 if (column >= 0 && column < 4)
1686                                         snprintf(buf, bufsize, "%s", sAngleParTitles[column]);
1687                                 return;
1688                         }
1689                         switch (column) {
1690                                 case 0: snprintf(buf, bufsize, "angle"); break;
1691                                 case 1:
1692                                         AtomTypeDecodeToString(ap->type1, types[0]);
1693                                         AtomTypeDecodeToString(ap->type2, types[1]);
1694                                         AtomTypeDecodeToString(ap->type3, types[2]);
1695                                         snprintf(buf, bufsize, "%s-%s-%s", types[0], types[1], types[2]);
1696                                         break;
1697                                 case 2:
1698                                 case 3:
1699                                         snprintf(buf, bufsize, "%.3f", (column == 2 ? ap->k * INTERNAL2KCAL : ap->a0 * kRad2Deg));
1700                                         break;
1701                         }
1702                         break;
1703                 }
1704                 case kDihedralParType: {
1705                         TorsionPar *tp = (TorsionPar *)up;
1706                         if (tp == NULL) {
1707                                 if (column >= 0 && column < 5)
1708                                         snprintf(buf, bufsize, "%s", sDihedralParTitles[column]);
1709                                 return;
1710                         }
1711                         switch (column) {
1712                                 case 0: snprintf(buf, bufsize, "dihe"); break;
1713                                 case 1:
1714                                         AtomTypeDecodeToString(tp->type1, types[0]);
1715                                         AtomTypeDecodeToString(tp->type2, types[1]);
1716                                         AtomTypeDecodeToString(tp->type3, types[2]);
1717                                         AtomTypeDecodeToString(tp->type4, types[3]);
1718                                         snprintf(buf, bufsize, "%s-%s-%s-%s", types[0], types[1], types[2], types[3]);
1719                                         break;
1720                                 case 3:
1721                                         snprintf(buf, bufsize, "%d", tp->period[0]);
1722                                         break;
1723                                 case 2:
1724                                 case 4:
1725                                         snprintf(buf, bufsize, "%.3f", (column == 2 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
1726                                         break;
1727                         }
1728                         break;
1729                 }
1730                 case kImproperParType: {
1731                         TorsionPar *tp = (TorsionPar *)up;
1732                         if (tp == NULL) {
1733                                 if (column >= 0 && column < 5)
1734                                         snprintf(buf, bufsize, "%s", sImproperParTitles[column]);
1735                                 return;
1736                         }
1737                         switch (column) {
1738                                 case 0: snprintf(buf, bufsize, "impr"); break;
1739                                 case 1:
1740                                         AtomTypeDecodeToString(tp->type1, types[0]);
1741                                         AtomTypeDecodeToString(tp->type2, types[1]);
1742                                         AtomTypeDecodeToString(tp->type3, types[2]);
1743                                         AtomTypeDecodeToString(tp->type4, types[3]);
1744                                         snprintf(buf, bufsize, "%s-%s-%s-%s", types[0], types[1], types[2], types[3]);
1745                                         break;
1746                                 case 3:
1747                                         snprintf(buf, bufsize, "%d", tp->period[0]);
1748                                         break;
1749                                 case 2:
1750                                 case 4:
1751                                         snprintf(buf, bufsize, "%.3f", (column == 2 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
1752                                         break;
1753                         }
1754                         break;
1755                 }
1756                 case kVdwPairParType: {
1757                         VdwPairPar *vp = (VdwPairPar *)up;
1758                         if (vp == NULL) {
1759                                 if (column >= 0 && column < 6)
1760                                         snprintf(buf, bufsize, "%s", sVdwPairParTitles[column]);
1761                                 return;
1762                         }
1763                         switch (column) {
1764                                 case 0: snprintf(buf, bufsize, "pvdw"); break;
1765                                 case 1:
1766                                         AtomTypeDecodeToString(vp->type1, types[0]);
1767                                         AtomTypeDecodeToString(vp->type2, types[1]);
1768                                         snprintf(buf, bufsize, "%s-%s", types[0], types[1]);
1769                                         break;
1770                                 case 2:
1771                                         snprintf(buf, bufsize, "%.6f", vp->eps * INTERNAL2KCAL);
1772                                         break;
1773                                 case 3:
1774                                         snprintf(buf, bufsize, "%.6f", vp->r_eq);
1775                                         break;
1776                                 case 4:
1777                                         snprintf(buf, bufsize, "%.6f", (vp->A14 == 0.0 ? 0.0 : vp->B14 * vp->B14 / vp->A14 * 0.25 * INTERNAL2KCAL));
1778                                         break;
1779                                 case 5:
1780                                         snprintf(buf, bufsize, "%.6f", vp->eps14 * INTERNAL2KCAL);
1781                                         break;
1782                         }
1783                         break;
1784                 }
1785 /*              case kElementParType: {
1786                         ElementPar *ap = (ElementPar *)up;
1787                         if (ap == NULL) {
1788                                 if (column >= 0 && column < 8)
1789                                         snprintf(buf, bufsize, "%s", sAtomParTitles[column]);
1790                                 return;
1791                         }
1792                         switch (column) {
1793                                 case 0: snprintf(buf, bufsize, "disp"); break;
1794                                 case 1: snprintf(buf, bufsize, "%.4s", ap->name); break;
1795                                 case 2: snprintf(buf, bufsize, "%d", ap->number); break;
1796                                 case 3: snprintf(buf, bufsize, "%.2f", ap->radius); break;
1797                                 case 4: snprintf(buf, bufsize, "%.3f", ap->r); break;
1798                                 case 5: snprintf(buf, bufsize, "%.3f", ap->g); break;
1799                                 case 6: snprintf(buf, bufsize, "%.3f", ap->b); break;
1800                                 case 7: snprintf(buf, bufsize, "%.3f", ap->weight); break;
1801                         }
1802                         break;
1803                 } */
1804                 default: return;
1805         }
1806         if (up != NULL && (column == 8 || column == 9)) {
1807                 if (column == 8 && ((BondPar *)up)->src == -1)
1808                         snprintf(buf, bufsize, "!NONE!");
1809                 else if ((p = ParameterGetComment(column == 8 ? ((BondPar *)up)->src : ((BondPar *)up)->com)) != NULL)
1810                         snprintf(buf, bufsize, "%s", p);
1811         }
1812 }
1813
1814 /*  Return values:
1815     -3: invalid
1816     -2: separator item
1817     -1: missing parameter
1818      0: molecule-local parameter
1819      1 and larger: global parameter values (gGlobalParInfo index + 1)  */
1820 int
1821 ParameterTableGetItemSource(Parameter *par, int row)
1822 {
1823         int src, type;
1824         UnionPar *up;
1825         if (par == NULL || row < 0)
1826                 return -3;
1827         up = ParameterTableGetItemPtr(par, row, &type);
1828         src = (type == kInvalidParType ? -3 : (up == NULL ? -2 : ((BondPar *)up)->src));
1829         if (type == kInvalidParType)
1830                 src = -3;
1831         else if (up == NULL)
1832                 src = -2;
1833         else src = ((BondPar *)up)->src;
1834         if (src > 0) {
1835                 /*  Search src in gGlobalParInfo  */
1836                 int i;
1837                 for (i = 0; i < gGlobalParInfo.count; i++) {
1838                         if (gGlobalParInfo.files[i].src == src)
1839                                 return i + 1;
1840                 }
1841                 return -3;  /*  Must not happen  */
1842         }
1843         return src;
1844 }
1845
1846 int
1847 ParameterTableIsItemEditable(Parameter *par, int column, int row)
1848 {
1849         UnionPar *up;
1850         int type, f;
1851         if (par == NULL || row < 0)
1852                 return 0;
1853         up = ParameterTableGetItemPtr(par, row, &type);
1854         if (type != kInvalidParType && up != NULL) {
1855                 /*  Valid type, not separator row, molecule-local value  */
1856                 int src = ((BondPar *)up)->src;
1857                 f = (src == 0 || src == -1);
1858         } else f = 0;
1859         switch (type) {
1860                 case kVdwParType: return (f && column > 0 && column < 8);
1861                 case kBondParType: return (f && column > 0 && column < 4);
1862                 case kAngleParType: return (f && column > 0 && column < 4);
1863                 case kDihedralParType: return (f && column > 0 && column < 5);
1864                 case kImproperParType: return (f && column > 0 && column < 5);
1865                 case kVdwPairParType: return (f && column > 0 && column < 5);
1866         /*      case kElementParType: return (f && column > 0 && column < 7); */
1867                 default: return 0;
1868         }
1869 }
1870
1871 #pragma mark ====== Utility Functions ======
1872
1873 int
1874 ElementParameterInitialize(const char *fname, char **outWarningMessage)
1875 {
1876         char buf[1024], name[6], fullname[16];
1877         float val[6];
1878         FILE *fp = NULL;
1879         int i, lineNumber, retval = 0;
1880         char *wbuf = NULL;
1881
1882         fp = fopen(fname, "rb");
1883         if (fp == NULL) {
1884                 retval = -1;
1885                 goto exit;
1886         }
1887         lineNumber = 0;
1888         while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
1889                 ElementPar *ep;
1890                 if (strncmp(buf, "element ", 8) != 0)
1891                         continue;  /*  Skip non-relevant lines  */
1892                 fullname[0] = 0;
1893                 if (sscanf(buf + 8, " %4s %f %f %f %f %f %f %15s", name, &val[0], &val[1], &val[2], &val[3], &val[4], &val[5], fullname) < 7) {
1894                         asprintf(&wbuf, "%s:%d: missing parameter in ELEMENT record", fname, lineNumber);
1895                         retval = 1;
1896                         goto exit;
1897                 }
1898                 i = (int)val[0];
1899                 if (i < 0 || i >= 200) {
1900                         asprintf(&wbuf, "%s:%d: The atomic number (%d) in ELEMENT record is out of range", fname, lineNumber, i);
1901                         retval = 2;
1902                         goto exit;
1903                 }
1904                 ep = AssignArray(&gElementParameters, &gCountElementParameters, sizeof(ElementPar), i, NULL);
1905                 memmove(ep->name, name, 4);
1906                 ep->number = i;
1907                 ep->radius = val[1];
1908                 ep->r = val[2];
1909                 ep->g = val[3];
1910                 ep->b = val[4];
1911                 ep->weight = val[5];
1912                 fullname[15] = 0;
1913                 memmove(ep->fullname, fullname, 16);
1914         }
1915 exit:
1916         if (fp != NULL)
1917                 fclose(fp);
1918         if (outWarningMessage != NULL)
1919                 *outWarningMessage = wbuf;
1920         return retval;
1921 }
1922
1923 UInt
1924 AtomTypeEncodeToUInt(const char *s)
1925 {
1926         Int i;
1927         UInt n, t;
1928         if ((s[0] == 'x' || s[0] == 'X') && s[1] == 0)
1929                 return kAtomTypeWildcard;
1930         if (s[0] >= '0' && s[0] <= '9')
1931                 return atoi(s);
1932         for (i = t = 0; i < 4; i++, s++) {
1933                 /*  Encode: variant*96*96*96*96 + a[0]*96*96*96 + a[1]*96*96 + a[2] * 96 + a[3]  */
1934                 static const UInt s_coeff[4] = {96*96*96, 96*96, 96, 1};
1935                 while (*s == ' ')
1936                         s++;
1937                 if (*s == 0)
1938                         break;
1939                 if (*s == ',' || *s == '.' || *s == ';' || *s == ':') {
1940                         /*  Variant (only [0-9a-z] are allowed) */
1941                         s++;
1942                         if (*s >= '0' && *s <= '9')
1943                                 n = *s - '0' + 1;
1944                         else if (*s >= 'A' && *s <= 'Z')
1945                                 n = *s - 'A' + 11;
1946                         else if (*s >= 'a' && *s <= 'z')
1947                                 n = *s - 'a' + 11;
1948                         else n = (*s % 36) + 1;  /*  Map to something non-zero  */
1949                         t += n * (96*96*96*96);
1950                         break;
1951                 }
1952                 n = (*s - 0x20) % 96;  /*  Map to 1..95  */
1953                 if (i == 0 && n < 32)
1954                         n = 32;
1955                 t += n * s_coeff[i];
1956         }
1957         return t;
1958 }
1959
1960 char *
1961 AtomTypeDecodeToString(UInt type, char *s)
1962 {
1963         static const UInt s_coeff[4] = {96*96*96, 96*96, 96, 1};
1964         static char buf[8];
1965         int i;
1966         UInt variant, n;
1967         if (s == NULL) {
1968                 s = buf;
1969                 buf[6] = 0;
1970         }
1971         if (type == kAtomTypeWildcard) {
1972                 s[0] = 'X';
1973                 s[1] = 0;
1974                 return s;
1975         }
1976         if (type < kAtomTypeMinimum) {
1977                 snprintf(s, 6, "%d", type);
1978                 return s;
1979         }
1980         for (i = 0; i < 4; i++) {
1981                 s[i] = (type / s_coeff[i]) % 96;
1982                 if (s[i] != 0)
1983                         s[i] += 0x20;
1984         }
1985         s[4] = 0;
1986         n = strlen(s);
1987         if ((variant = (type / (96*96*96*96)) % 96) != 0) {
1988                 s[n] = '.';
1989                 s[n + 1] = (variant <= 10 ? '0' + variant - 1 : 'a' + variant - 11);
1990                 s[n + 2] = 0;
1991         }
1992         return s;
1993 }
1994
1995 Int
1996 ElementToInt(const char *s)
1997 {
1998         int i;
1999         ElementPar *p;
2000         for (i = 0, p = gElementParameters; i < gCountElementParameters; i++, p++) {
2001                 if (p->name[0] == toupper(s[0]) && p->name[1] == tolower(s[1]))
2002                         return p->number;
2003         }
2004         return 0;
2005 }
2006
2007 char *
2008 ElementToString(Int elem, char *s)
2009 {
2010         if (elem >= 0 && elem < gCountElementParameters) {
2011                 const char *cs = gElementParameters[elem].name;
2012                 s[0] = cs[0];
2013                 s[1] = cs[1];
2014                 s[2] = 0;
2015                 return s;
2016         } else return NULL;
2017 }
2018
2019 Int
2020 AtomNameToElement(const char *s)
2021 {
2022         char element[4];
2023         const char *p;
2024         /*  $element = ($name =~ /([A-Za-z]{1,2})/); # in Perl  */
2025         element[0] = 0;
2026         for (p = s; *p != 0; p++) {
2027                 if (isalpha(*p) && *p != '_') {
2028                         element[0] = toupper(*p);
2029                         if (isalpha(p[1]) && p[1] != '_') {
2030                                 element[1] = toupper(p[1]);
2031                                 element[2] = 0;
2032                         } else {
2033                                 element[1] = 0;
2034                         }
2035                         break;
2036                 }
2037         }
2038         return ElementToInt(element);
2039 }
2040
2041 Int
2042 GuessAtomicNumber(const char *name, Double weight)
2043 {
2044         int i;
2045         ElementPar *p;
2046         char buf[4];
2047         const char *cp;
2048         for (i = 0, p = gElementParameters; i < gCountElementParameters; i++, p++) {
2049                 if (p->weight > 0.0 && fabs(weight - p->weight) < 0.1)
2050                         return p->number;
2051         }
2052         buf[0] = 0;
2053         for (cp = name; *cp != 0 && cp < name + 4; cp++) {
2054                 if (isalpha(*cp) && *cp != '_') {
2055                         buf[0] = toupper(*cp);
2056                         if (isalpha(cp[1]) && cp[1] != '_') {
2057                                 buf[1] = toupper(cp[1]);
2058                                 buf[2] = 0;
2059                         } else {
2060                                 buf[1] = 0;
2061                         }
2062                         break;
2063                 }
2064         }
2065         return ElementToInt(buf);
2066 }
2067
2068 Double
2069 WeightForAtomicNumber(Int elem)
2070 {
2071         if (elem >= 1 && elem < gCountElementParameters)
2072                 return gElementParameters[elem].weight;
2073         else return 0.0;
2074 }
2075
2076 Double
2077 RadiusForAtomicNumber(Int elem)
2078 {
2079         if (elem >= 1 && elem < gCountElementParameters)
2080                 return gElementParameters[elem].radius;
2081         else return 0.0;
2082 }
2083