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