4 * Created by Toshi Nagata on 06/03/11.
5 * Copyright 2006-2008 Toshi Nagata. All rights reserved.
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.
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.
23 /* Global parameter: it is initialized by the first call to ParameterReadFromFile() */
24 Parameter *gBuiltinParameters = NULL;
26 /* Global parameter */
27 AtomPar *gDispAtomParameters = NULL;
28 Int gCountDispAtomParameters = 0;
30 static Parameter *sParameterRoot = NULL;
31 static int sParameterUntitledCount = 0;
33 /* Global struct for storing information for global parameter files */
34 GlobalParInfoRecord gGlobalParInfo;
36 #pragma mark ====== Parameter Alloc/Release ======
42 Parameter *par = (Parameter *)calloc(sizeof(Parameter), 1);
44 Panic("Cannot allocate new parameter record");
45 snprintf(name, sizeof name, "Untitled %d", sParameterUntitledCount++);
46 ObjectInit((Object *)par, (Object **)&sParameterRoot, name);
51 ParameterDuplicate(const Parameter *par)
53 Parameter *npar = ParameterNew();
54 if (par->bondPars != NULL) {
55 npar->bondPars = (BondPar *)malloc(sizeof(BondPar) * par->nbondPars);
56 if (npar->bondPars == NULL)
58 memmove(npar->bondPars, par->bondPars, sizeof(BondPar) * par->nbondPars);
59 npar->nbondPars = par->nbondPars;
61 if (par->anglePars != NULL) {
62 npar->anglePars = (AnglePar *)malloc(sizeof(AnglePar) * par->nanglePars);
63 if (npar->anglePars == NULL)
65 memmove(npar->anglePars, par->anglePars, sizeof(AnglePar) * par->nanglePars);
66 npar->nanglePars = par->nanglePars;
68 if (par->dihedralPars != NULL) {
69 npar->dihedralPars = (TorsionPar *)malloc(sizeof(TorsionPar) * par->ndihedralPars);
70 if (npar->dihedralPars == NULL)
72 memmove(npar->dihedralPars, par->dihedralPars, sizeof(TorsionPar) * par->ndihedralPars);
73 npar->ndihedralPars = par->ndihedralPars;
75 if (par->improperPars != NULL) {
76 npar->improperPars = (TorsionPar *)malloc(sizeof(TorsionPar) * par->nimproperPars);
77 if (npar->improperPars == NULL)
79 memmove(npar->improperPars, par->improperPars, sizeof(TorsionPar) * par->nimproperPars);
80 npar->nimproperPars = par->nimproperPars;
82 if (par->vdwPars != NULL) {
83 npar->vdwPars = (VdwPar *)malloc(sizeof(VdwPar) * par->nvdwPars);
84 if (npar->vdwPars == NULL)
86 memmove(npar->vdwPars, par->vdwPars, sizeof(VdwPar) * par->nvdwPars);
87 npar->nvdwPars = par->nvdwPars;
89 if (par->vdwpPars != NULL) {
90 npar->vdwpPars = (VdwPairPar *)malloc(sizeof(VdwPairPar) * par->nvdwpPars);
91 if (npar->vdwpPars == NULL)
93 memmove(npar->vdwpPars, par->vdwpPars, sizeof(VdwPairPar) * par->nvdwpPars);
94 npar->nvdwpPars = par->nvdwpPars;
96 if (par->vdwCutoffPars != NULL) {
97 npar->vdwCutoffPars = (VdwCutoffPar *)malloc(sizeof(VdwCutoffPar) * par->nvdwCutoffPars);
98 if (npar->vdwCutoffPars == NULL)
100 memmove(npar->vdwCutoffPars, par->vdwCutoffPars, sizeof(VdwCutoffPar) * par->nvdwCutoffPars);
101 npar->nvdwCutoffPars = par->nvdwCutoffPars;
103 /* if (par->atomPars != NULL) {
104 npar->atomPars = (AtomPar *)malloc(sizeof(AtomPar) * par->natomPars);
105 if (npar->atomPars == NULL)
107 memmove(npar->atomPars, par->atomPars, sizeof(AtomPar) * par->natomPars);
108 npar->natomPars = par->natomPars;
112 ParameterRelease(npar);
117 ParameterWithName(const char *name)
119 return (Parameter *)ObjectWithName(name, (Object *)sParameterRoot);
122 /* Assign a unique name to this parameter record */
124 ParameterSetName(Parameter *par, const char *name)
126 ObjectSetName((Object *)par, name, (Object *)sParameterRoot);
130 ParameterGetName(Parameter *par)
132 return ObjectGetName((Object *)par);
136 ParameterRetain(Parameter *par)
138 ObjectIncrRefCount((Object *)par);
142 ParameterRelease(Parameter *par)
144 if (ObjectDecrRefCount((Object *)par) == 0) {
145 if (par->bondPars != NULL)
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)
155 if (par->vdwpPars != NULL)
157 if (par->vdwCutoffPars != NULL)
158 free(par->vdwCutoffPars);
159 /* if (par->atomPars != NULL)
160 free(par->atomPars); */
161 ObjectDealloc((Object *)par, (Object **)&sParameterRoot);
165 #pragma mark ====== ParameterRef Definitions ======
168 ParameterRefNew(struct Molecule *mol, int type, int idx)
170 ParameterRef *pref = (ParameterRef *)calloc(sizeof(ParameterRef), 1);
175 pref->parType = type;
182 ParameterRefRelease(ParameterRef *pref)
185 if (pref->mol != NULL)
186 MoleculeRelease(pref->mol);
192 ParameterGetUnionParFromTypeAndIndex(Parameter *par, int type, int index)
198 if (index >= 0 && index < par->nbondPars)
199 return (UnionPar *)(par->bondPars + index);
202 if (index >= 0 && index < par->nanglePars)
203 return (UnionPar *)(par->anglePars + index);
205 case kDihedralParType:
206 if (index >= 0 && index < par->ndihedralPars)
207 return (UnionPar *)(par->dihedralPars + index);
209 case kImproperParType:
210 if (index >= 0 && index < par->nimproperPars)
211 return (UnionPar *)(par->improperPars + index);
214 if (index >= 0 && index < par->nvdwPars)
215 return (UnionPar *)(par->vdwPars + index);
217 case kVdwPairParType:
218 if (index >= 0 && index < par->nvdwpPars)
219 return (UnionPar *)(par->vdwpPars + index);
221 case kVdwCutoffParType:
222 if (index >= 0 && index < par->nvdwCutoffPars)
223 return (UnionPar *)(par->vdwCutoffPars + index);
225 /* case kAtomParType:
226 if (index >= 0 && index < par->natomPars)
227 return (UnionPar *)(par->atomPars + index);
235 ParameterGetCountForType(Parameter *par, int type)
241 return par->nbondPars;
243 return par->nanglePars;
244 case kDihedralParType:
245 return par->ndihedralPars;
246 case kImproperParType:
247 return par->nimproperPars;
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; */
262 ParameterGetSizeForType(int type)
266 return sizeof(BondPar);
268 return sizeof(AnglePar);
269 case kDihedralParType:
270 return sizeof(TorsionPar);
271 case kImproperParType:
272 return sizeof(TorsionPar);
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; */
287 ParameterRefGetPar(ParameterRef *pref)
292 if (pref->mol != NULL)
293 par = pref->mol->par;
294 else par = gBuiltinParameters;
297 return ParameterGetUnionParFromTypeAndIndex(par, pref->parType, pref->idx);
300 #pragma mark ====== Insert/Delete (for MolAction) ======
303 ParameterInsert(Parameter *par, Int type, const UnionPar *up, struct IntGroup *where)
305 Int i, n1, n2, size, *ip;
312 ip = &par->nbondPars;
313 size = sizeof(BondPar);
317 ip = &par->nanglePars;
318 size = sizeof(AnglePar);
320 case kDihedralParType:
321 p = &par->dihedralPars;
322 ip = &par->ndihedralPars;
323 size = sizeof(TorsionPar);
325 case kImproperParType:
326 p = &par->improperPars;
327 ip = &par->nimproperPars;
328 size = sizeof(TorsionPar);
333 size = sizeof(VdwPar);
335 case kVdwPairParType:
337 ip = &par->nvdwpPars;
338 size = sizeof(VdwPairPar);
340 case kVdwCutoffParType:
341 p = &par->vdwCutoffPars;
342 ip = &par->nvdwCutoffPars;
343 size = sizeof(VdwCutoffPar);
345 /* case kAtomParType:
347 ip = &par->natomPars;
348 size = sizeof(AtomPar);
354 for (i = 0; (n1 = IntGroupGetNthPoint(where, i)) >= 0; i++) {
355 if (InsertArray(p, ip, size, n1, 1, up + i) == NULL)
363 sParameterDeleteOrCopy(Parameter *par, Int type, UnionPar *up, struct IntGroup *where, Int copyflag)
365 Int i, n1, n2, size, *ip;
371 p = (char **)&par->bondPars;
372 ip = &par->nbondPars;
373 size = sizeof(BondPar);
376 p = (char **)&par->anglePars;
377 ip = &par->nanglePars;
378 size = sizeof(AnglePar);
380 case kDihedralParType:
381 p = (char **)&par->dihedralPars;
382 ip = &par->ndihedralPars;
383 size = sizeof(TorsionPar);
385 case kImproperParType:
386 p = (char **)&par->improperPars;
387 ip = &par->nimproperPars;
388 size = sizeof(TorsionPar);
391 p = (char **)&par->vdwPars;
393 size = sizeof(VdwPar);
395 case kVdwPairParType:
396 p = (char **)&par->vdwpPars;
397 ip = &par->nvdwpPars;
398 size = sizeof(VdwPairPar);
400 case kVdwCutoffParType:
401 p = (char **)&par->vdwCutoffPars;
402 ip = &par->nvdwCutoffPars;
403 size = sizeof(VdwCutoffPar);
405 /* case kAtomParType:
406 p = (char **)&par->atomPars;
407 ip = &par->natomPars;
408 size = sizeof(AtomPar);
414 for (i = IntGroupGetCount(where) - 1; i >= 0 && (n1 = IntGroupGetNthPoint(where, i)) >= 0; i--) {
416 return -1; /* Internal error */
419 memmove(up + i, *p + size * n1, size);
421 if (DeleteArray(p, ip, size, n1, 1, (up ? up + i : NULL)) == NULL)
430 ParameterDelete(Parameter *par, Int type, UnionPar *up, struct IntGroup *where)
432 return sParameterDeleteOrCopy(par, type, up, where, 0);
436 ParameterCopy(Parameter *par, Int type, UnionPar *up, struct IntGroup *where)
438 return sParameterDeleteOrCopy(par, type, up, where, 1);
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. */
444 ParameterRenumberAtoms(Int type, UnionPar *up, Int oldnatoms, const Int *old2new)
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)))
453 return RENUMBER_FIELD(BondPar, type1) + RENUMBER_FIELD(BondPar, type2);
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);
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);
472 #pragma mark ====== Load from Files ======
475 s_AppendWarning(char **ptr, const char *fmt, ...)
483 *ptr = (char *)malloc(128);
491 for (len = 128; len <= n; len *= 2);
494 nn = vasprintf(&buf, fmt, ap);
498 while (len <= n + nn)
500 *ptr = (char *)realloc(*ptr, len);
504 strncpy(*ptr + n, buf, nn);
505 *(*ptr + n + nn) = 0;
511 ParameterGlobalParIndexForSrcIndex(int src)
514 if (src == 0 || src == -1)
516 for (i = 0; i < gGlobalParInfo.count; i++) {
517 if (gGlobalParInfo.files[i].src == src)
524 ParameterCommentIndexForGlobalFileName(const char *p)
526 const char *pp = p + strlen(p), *p1 = NULL, *p2 = NULL;
529 if (p2 == NULL && *pp == '.')
531 if (p1 == NULL && (*pp == '\'' || *pp == '/'))
538 snprintf(buf, sizeof(buf), "%.*s", (int)(p2 - p1), p1);
539 return ParameterCommentIndex(buf);
543 ParameterCompare(const UnionPar *up1, const UnionPar *up2, int type)
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);
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);
561 case kDihedralParType:
562 case kImproperParType: {
563 const TorsionPar *tp1 = &up1->torsion;
564 const TorsionPar *tp2 = &up2->torsion;
566 if (tp1->mult != tp2->mult)
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))
573 if (tp1->type3 != tp2->type3)
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))
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)
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);
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);
599 case kVdwCutoffParType: {
600 const VdwCutoffPar *vp1 = &up1->vdwcutoff;
601 const VdwCutoffPar *vp2 = &up2->vdwcutoff;
602 if (vp1->type != vp2->type)
604 if (vp1->type == 0) {
605 if ((vp1->n1 != vp2->n1 || vp1->n2 != vp2->n2) && (vp1->n1 != vp2->n2 || vp1->n2 != vp2->n1))
608 if (vp1->n1 != vp2->n1 || vp1->n2 != vp2->n2 || vp1->n3 != vp2->n3 || vp1->n4 != vp2->n4)
611 return fabs(vp1->cutoff - vp2->cutoff) < 1e-8;
617 /* bp can also be AnglePar *, TorsionPar *, etc. */
619 s_StoreComment(int parType, BondPar *bp, char *p, const char *src)
621 char *s, *pp, *pcom, buf[40];
622 int embedded = 0, src_idx;
625 while (isspace(*p) || *p == ';' || *p == '!')
628 if (src == NULL && *p == '[') {
629 /* Source is embedded */
633 while (*p != ']' && *p != 0)
636 if (len >= sizeof(buf))
637 len = sizeof(buf) - 1;
638 strncpy(buf, s, len);
649 while (*pp != 0 && *pp != '\r' && *pp != '\n')
652 if (src != NULL && *src != 0) {
653 src_idx = ParameterCommentIndex(src);
655 /* Compare with already known global parameters, and src is set if the same one is found */
658 for (i = 0; (up2 = ParameterGetUnionParFromTypeAndIndex(gBuiltinParameters, parType, i)) != NULL; i++) {
659 if (up2->bond.src == src_idx && ParameterCompare((UnionPar *)bp, up2, parType)) {
665 /* Not found: embedded source is retained, and this entry is regarded as "no source" */
669 } else bp->src = src_idx;
672 bp->com = ParameterCommentIndex(p);
675 /* bp can also be BondPar*, AnglePar *, TorsionPar *, etc. */
677 s_CommentToString(char *buf, int bufsize, void *bp)
679 const char *src, *com;
680 src = ParameterGetComment(((BondPar *)bp)->src);
681 com = ParameterGetComment(((BondPar *)bp)->com);
682 if (src == NULL && com == NULL)
684 else if (src != NULL)
685 snprintf(buf, bufsize, " ![%s] %s", src, (com == NULL ? "" : com));
687 snprintf(buf, bufsize, " ! %s", com);
691 ParameterReadFromString(Parameter *par, char *buf, char **wbufp, const char *fname, int lineNumber, int src_idx)
694 char com[12], com2[12], type[4][8];
696 float val[6]; /* not Double */
701 if (sscanf(buf, " %11s", com) <= 0 || !isalpha(com[0]))
704 for (i = 0; i < len; i++)
705 com[i] = tolower(com[i]);
706 if (strncmp(com, "include", len) == 0) {
711 for ( ; (c = buf[len]) != 0; len++) {
712 if (c == ' ' || c == '\t')
730 p = (char *)malloc(strlen(fname) + i + 1);
735 pp = strrchr(p, '/');
740 i = ParameterReadFromFile(par, p, &wbuf1, &wc1);
742 s_AppendWarning(wbufp, "In included file %s:\n%s", p, wbuf1);
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. */
759 else src = ParameterGetComment(src_idx);
760 options = kParameterLookupNoWildcard | kParameterLookupNoBaseAtomType;
761 if (strncmp(com, "bond", len) == 0) {
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);
767 itype[0] = AtomTypeEncodeToUInt(type[0]);
768 itype[1] = AtomTypeEncodeToUInt(type[1]);
769 if (itype[0] > itype[1]) {
774 val[0] *= KCAL2INTERNAL;
775 bp = ParameterLookupBondPar(par, itype[0], itype[1], options);
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]);
782 bp = AssignArray(&par->bondPars, &par->nbondPars, sizeof(*bp), par->nbondPars, NULL);
783 bp->type1 = itype[0];
784 bp->type2 = itype[1];
787 s_StoreComment(kBondParType, bp, buf + n, src);
788 } else if (strncmp(com, "angle", len) == 0) {
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);
794 itype[0] = AtomTypeEncodeToUInt(type[0]);
795 itype[1] = AtomTypeEncodeToUInt(type[1]);
796 itype[2] = AtomTypeEncodeToUInt(type[2]);
797 if (itype[0] > itype[2]) {
802 val[0] *= KCAL2INTERNAL;
803 val[1] *= (3.14159265358979 / 180.0);
804 ap = ParameterLookupAnglePar(par, itype[0], itype[1], itype[2], options);
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]);
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];
817 s_StoreComment(kAngleParType, (BondPar *)ap, buf + n, src);
818 } else if (strncmp(com, "dihedral", len) == 0) {
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);
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]) {
836 dp = ParameterLookupDihedralPar(par, itype[0], itype[1], itype[2], itype[3], options);
837 val[0] *= KCAL2INTERNAL;
838 val[1] *= 3.14159265358979 / 180.0;
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]);
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];
851 dp->period[0] = ival[0];
852 dp->phi0[0] = val[1];
854 s_StoreComment(kDihedralParType, (BondPar *)dp, buf + n, src);
855 } else if (strncmp(com, "improper", len) == 0) {
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);
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]) {
870 if (itype[0] > itype[3]) {
875 if (itype[1] > itype[3]) {
880 ip = ParameterLookupImproperPar(par, itype[0], itype[1], itype[2], itype[3], options);
881 val[0] *= KCAL2INTERNAL;
882 val[1] *= 3.14159265358979 / 180.0;
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]);
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];
895 ip->period[0] = ival[0];
896 ip->phi0[0] = val[1];
898 s_StoreComment(kImproperParType, (BondPar *)ip, buf + n, src);
899 } else if (strncmp(com, "nonbonded", len) == 0 || strncmp(com, "vdw", len) == 0) {
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"));
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;
922 val[4] = val[5] = 0.0;
923 if (sscanf(p, "%d %n", &ival[0], &n) == 1) {
924 vtemp.atomicNumber = ival[0];
927 if (sscanf(p, "%f %n", &val[4], &n) == 1) {
928 vtemp.weight = val[4];
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];
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;
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]);
957 vp = AssignArray(&par->vdwPars, &par->nvdwPars, sizeof(*vp), i, NULL);
958 vtemp.com = vp->com; /* Keep comment field */
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"));
968 itype[0] = AtomTypeEncodeToUInt(type[0]);
969 itype[1] = AtomTypeEncodeToUInt(type[1]);
970 if (itype[0] > itype[1]) {
975 vtemp.type1 = itype[0];
976 vtemp.type2 = itype[1];
977 if (flag) { /* eps/r_eq representation */
978 vtemp.eps = val[0] * KCAL2INTERNAL;
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;
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]);
1007 vp = AssignArray(&par->vdwpPars, &par->nvdwpPars, sizeof(*vp), i, NULL);
1008 vtemp.com = vp->com; /* Keep comment field */
1010 s_StoreComment(kVdwPairParType, (BondPar *)vp, buf + n, src);
1012 s_AppendWarning(wbufp, "%s:%d: unknown keyword %s\n", fname, lineNumber, com);
1019 ParameterReadFromFile(Parameter *par, const char *fname, char **outWarningMessage, int *outWarningCount)
1030 par = gBuiltinParameters;
1032 par = ParameterNew();
1033 gBuiltinParameters = par;
1041 fp = fopen(fname, "rb");
1043 s_AppendWarning(&wbuf, "Cannot open parameter file %s\n", fname);
1048 if (par != gBuiltinParameters || first)
1051 src_idx = ParameterCommentIndexForGlobalFileName(fname);
1053 if (par == gBuiltinParameters && !first) {
1054 /* Ensure the "source" field is unique */
1058 for (i = 0; i < gGlobalParInfo.count; i++) {
1059 if (gGlobalParInfo.files[i].src == src_idx)
1062 if (i < gGlobalParInfo.count) {
1063 /* Delete the existing Parameters from the same source */
1064 ParameterDeleteAllEntriesForSource(par, src_idx);
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 == '\\')
1077 ip->dir = (char *)malloc(i + 1);
1078 if (ip->dir != NULL) {
1079 strncpy(ip->dir, fname, i);
1084 ip->name = strdup(p);
1089 while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
1090 int wc1 = ParameterReadFromString(par, buf, &wbuf, fname, lineNumber, src_idx);
1099 gGlobalParInfo.builtinCount = gGlobalParInfo.count;
1104 if (outWarningMessage != NULL)
1105 *outWarningMessage = wbuf;
1106 else if (wbuf != NULL)
1108 if (outWarningCount != NULL)
1109 *outWarningCount = (wcount >= 0 ? wcount : 0);
1111 return (wcount >= 0 ? 0 : wcount);
1115 ParameterAppendToFile(Parameter *par, FILE *fp)
1120 int bufsize = sizeof(buf);
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);
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);
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);
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);
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 */
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);
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);
1208 ParameterDeleteAllEntriesForSource(Parameter *par, int src_idx)
1213 /* if (fname == NULL)
1215 src_idx = ParameterCommentIndexForGlobalFileName(fname); */
1219 for (type = kFirstParType; type <= kLastParType; type++) {
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);
1226 i = IntGroupGetCount(ig);
1228 ParameterDelete(par, type, NULL, ig);
1231 IntGroupRelease(ig);
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);
1245 #pragma mark ====== Parameter Comments ======
1247 static const char **sParameterComments;
1248 static Int sNumberOfParameterComments;
1251 ParameterCommentIndex(const char *comment)
1255 if (comment == NULL || comment[0] == 0)
1257 /* Duplicate the comment, convert whitespaces to ' ', and chop trailing whitespaces */
1258 p = strdup(comment);
1260 for (pp = p + strlen(p) - 1; pp >= p; pp--) {
1267 for (i = 1; i < sNumberOfParameterComments; i++) {
1268 if (strcmp(sParameterComments[i], p) == 0) {
1273 if (sNumberOfParameterComments == 0) {
1274 /* Index 0 is skipped */
1275 AssignArray(&sParameterComments, &sNumberOfParameterComments, sizeof(char *), 1, &p);
1277 AssignArray(&sParameterComments, &sNumberOfParameterComments, sizeof(char *), sNumberOfParameterComments, &p);
1279 return sNumberOfParameterComments - 1;
1283 ParameterGetComment(int idx)
1285 if (idx <= 0 || idx >= sNumberOfParameterComments)
1286 return NULL; /* No such number */
1287 return sParameterComments[idx];
1290 #pragma mark ====== Parameter Lookup ======
1292 #define s_ParMatch(_t1, _t2, _nowildcard) (_t1 == _t2 || (!_nowildcard && _t1 == kAtomTypeWildcard))
1294 /* Returns non-zero if the parameter record contains designated atom_type.
1295 The atom_type can also be an atom index. */
1297 ParameterDoesContainAtom(Int type, UnionPar *up, UInt atom_type, Int options)
1299 #define CHECK_FIELD(_tp, _field) s_ParMatch((((_tp *)up)->_field), atom_type, nowildcard)
1300 Int nowildcard = (options & kParameterLookupNoWildcard);
1303 return CHECK_FIELD(BondPar, type1) || CHECK_FIELD(BondPar, type2);
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);
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);
1322 ParameterLookupBondPar(Parameter *par, UInt t1, UInt t2, Int options)
1326 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
1336 if (s_ParMatch(bp->type1, t1, nowildcard) && s_ParMatch(bp->type2, t2, nowildcard))
1338 if (s_ParMatch(bp->type1, t2, nowildcard) && s_ParMatch(bp->type2, t1, nowildcard))
1341 if (options & kParameterLookupNoBaseAtomType)
1343 return ParameterLookupBondPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1347 ParameterLookupAnglePar(Parameter *par, UInt t1, UInt t2, UInt t3, Int options)
1351 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
1361 if (s_ParMatch(ap->type1, t1, nowildcard) && s_ParMatch(ap->type2, t2, nowildcard) && s_ParMatch(ap->type3, t3, nowildcard))
1363 if (s_ParMatch(ap->type1, t3, nowildcard) && s_ParMatch(ap->type2, t2, nowildcard) && s_ParMatch(ap->type3, t1, nowildcard))
1366 if (options & kParameterLookupNoBaseAtomType)
1368 return ParameterLookupAnglePar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1372 ParameterLookupDihedralPar(Parameter *par, UInt t1, UInt t2, UInt t3, UInt t4, Int options)
1376 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
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))
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))
1391 if (options & kParameterLookupNoBaseAtomType)
1393 return ParameterLookupDihedralPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, t4 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1397 ParameterLookupImproperPar(Parameter *par, UInt t1, UInt t2, UInt t3, UInt t4, Int options)
1401 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
1411 if (!s_ParMatch(tp->type3, t3, nowildcard))
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)))
1421 if (options & kParameterLookupNoBaseAtomType)
1423 return ParameterLookupImproperPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, t3 % kAtomTypeVariantBase, t4 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1427 ParameterLookupVdwPar(Parameter *par, UInt t1, Int options)
1431 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
1441 if (s_ParMatch(vp->type1, t1, nowildcard))
1444 if (options & kParameterLookupNoBaseAtomType)
1446 return ParameterLookupVdwPar(par, t1 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1450 ParameterLookupVdwPairPar(Parameter *par, UInt t1, UInt t2, Int options)
1454 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
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)))
1468 if (options & kParameterLookupNoBaseAtomType)
1470 return ParameterLookupVdwPairPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1474 ParameterLookupVdwCutoffPar(Parameter *par, UInt t1, UInt t2, Int options)
1478 Int nowildcard = (options & kParameterLookupNoWildcard);
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)))
1488 if (vp->type == 0) {
1489 if (s_ParMatch(vp->n1, t1, nowildcard) && s_ParMatch(vp->n2, t2, nowildcard))
1491 if (s_ParMatch(vp->n1, t2, nowildcard) && s_ParMatch(vp->n2, t1, nowildcard))
1494 if (vp->n1 <= t1 && vp->n2 >= t1 && vp->n3 <= t2 && vp->n4 <= t2)
1496 if (vp->n1 <= t2 && vp->n2 >= t2 && vp->n3 <= t1 && vp->n4 >= t1)
1500 if (options & kParameterLookupNoBaseAtomType)
1502 return ParameterLookupVdwCutoffPar(par, t1 % kAtomTypeVariantBase, t2 % kAtomTypeVariantBase, options | kParameterLookupNoBaseAtomType);
1505 #pragma mark ====== Table View Support ======
1508 ParameterTableNumberOfRows(Parameter *par)
1512 return par->nvdwPars + par->nbondPars + par->nanglePars + par->ndihedralPars + par->nimproperPars + par->nvdwpPars + 6;
1516 ParameterTableGetItemIndex(Parameter *par, int row, int *type)
1518 if (par == NULL || row < 0) {
1519 *type = kInvalidParType;
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;
1535 *type = kInvalidParType;
1542 ParameterTableGetItemPtr(Parameter *par, int row, int *type)
1544 if (par == NULL || row < 0) {
1545 *type = kInvalidParType;
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); */
1570 *type = kInvalidParType;
1576 ParameterTableGetItemText(Parameter *par, int column, int row, char *buf, int bufsize)
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"}; */
1587 UnionPar *up = NULL;
1590 if (par == NULL || row < 0)
1592 up = ParameterTableGetItemPtr(par, row, &type);
1595 VdwPar *vp = (VdwPar *)up;
1597 if (column >= 0 && column < 8)
1598 snprintf(buf, bufsize, "%s", sVdwParTitles[column]);
1602 case 0: snprintf(buf, bufsize, "vdw"); break;
1604 AtomTypeDecodeToString(vp->type1, types[0]);
1605 snprintf(buf, bufsize, "%s", types[0]);
1608 snprintf(buf, bufsize, "%.5f", vp->eps * INTERNAL2KCAL);
1611 snprintf(buf, bufsize, "%.5f", vp->r_eq);
1614 snprintf(buf, bufsize, "%.5f", vp->eps14 * INTERNAL2KCAL);
1617 snprintf(buf, bufsize, "%.5f", vp->r_eq14);
1620 snprintf(buf, bufsize, "%d", vp->atomicNumber);
1623 snprintf(buf, bufsize, "%.3f", vp->weight);
1628 case kBondParType: {
1629 BondPar *bp = (BondPar *)up;
1631 if (column >= 0 && column < 4)
1632 snprintf(buf, bufsize, "%s", sBondParTitles[column]);
1636 case 0: snprintf(buf, bufsize, "bond"); break;
1638 AtomTypeDecodeToString(bp->type1, types[0]);
1639 AtomTypeDecodeToString(bp->type2, types[1]);
1640 snprintf(buf, bufsize, "%s-%s", types[0], types[1]);
1644 snprintf(buf, bufsize, "%.3f", (column == 2 ? bp->k * INTERNAL2KCAL : bp->r0));
1649 case kAngleParType: {
1650 AnglePar *ap = (AnglePar *)up;
1652 if (column >= 0 && column < 4)
1653 snprintf(buf, bufsize, "%s", sAngleParTitles[column]);
1657 case 0: snprintf(buf, bufsize, "angle"); break;
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]);
1666 snprintf(buf, bufsize, "%.3f", (column == 2 ? ap->k * INTERNAL2KCAL : ap->a0 * kRad2Deg));
1671 case kDihedralParType: {
1672 TorsionPar *tp = (TorsionPar *)up;
1674 if (column >= 0 && column < 5)
1675 snprintf(buf, bufsize, "%s", sDihedralParTitles[column]);
1679 case 0: snprintf(buf, bufsize, "dihe"); break;
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]);
1688 snprintf(buf, bufsize, "%d", tp->period[0]);
1692 snprintf(buf, bufsize, "%.3f", (column == 2 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
1697 case kImproperParType: {
1698 TorsionPar *tp = (TorsionPar *)up;
1700 if (column >= 0 && column < 5)
1701 snprintf(buf, bufsize, "%s", sImproperParTitles[column]);
1705 case 0: snprintf(buf, bufsize, "impr"); break;
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]);
1714 snprintf(buf, bufsize, "%d", tp->period[0]);
1718 snprintf(buf, bufsize, "%.3f", (column == 2 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
1723 case kVdwPairParType: {
1724 VdwPairPar *vp = (VdwPairPar *)up;
1726 if (column >= 0 && column < 6)
1727 snprintf(buf, bufsize, "%s", sVdwPairParTitles[column]);
1731 case 0: snprintf(buf, bufsize, "pvdw"); break;
1733 AtomTypeDecodeToString(vp->type1, types[0]);
1734 AtomTypeDecodeToString(vp->type2, types[1]);
1735 snprintf(buf, bufsize, "%s-%s", types[0], types[1]);
1738 snprintf(buf, bufsize, "%.6f", vp->eps * INTERNAL2KCAL);
1741 snprintf(buf, bufsize, "%.6f", vp->r_eq);
1744 snprintf(buf, bufsize, "%.6f", (vp->A14 == 0.0 ? 0.0 : vp->B14 * vp->B14 / vp->A14 * 0.25 * INTERNAL2KCAL));
1747 snprintf(buf, bufsize, "%.6f", vp->eps14 * INTERNAL2KCAL);
1752 /* case kAtomParType: {
1753 AtomPar *ap = (AtomPar *)up;
1755 if (column >= 0 && column < 8)
1756 snprintf(buf, bufsize, "%s", sAtomParTitles[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;
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);
1784 -1: missing parameter
1785 0: molecule-local parameter
1786 1 and larger: global parameter values (gGlobalParInfo index + 1) */
1788 ParameterTableGetItemSource(Parameter *par, int row)
1792 if (par == NULL || row < 0)
1794 up = ParameterTableGetItemPtr(par, row, &type);
1795 src = (type == kInvalidParType ? -3 : (up == NULL ? -2 : ((BondPar *)up)->src));
1796 if (type == kInvalidParType)
1798 else if (up == NULL)
1800 else src = ((BondPar *)up)->src;
1802 /* Search src in gGlobalParInfo */
1804 for (i = 0; i < gGlobalParInfo.count; i++) {
1805 if (gGlobalParInfo.files[i].src == src)
1808 return -3; /* Must not happen */
1814 ParameterTableIsItemEditable(Parameter *par, int column, int row)
1818 if (par == NULL || row < 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);
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); */
1838 #pragma mark ====== Utility Functions ======
1841 AtomParameterInitialize(const char *fname, char **outWarningMessage)
1843 char buf[1024], name[6], fullname[16];
1846 int i, lineNumber, retval = 0;
1849 fp = fopen(fname, "rb");
1855 while (ReadLine(buf, sizeof buf, fp, &lineNumber) > 0) {
1857 if (strncmp(buf, "dispatom ", 9) != 0)
1858 continue; /* Skip non-relevant lines */
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);
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);
1871 ap = AssignArray(&gDispAtomParameters, &gCountDispAtomParameters, sizeof(AtomPar), i, NULL);
1872 memmove(ap->name, name, 4);
1874 ap->radius = val[1];
1878 ap->weight = val[5];
1880 memmove(ap->fullname, fullname, 16);
1885 if (outWarningMessage != NULL)
1886 *outWarningMessage = wbuf;
1891 AtomTypeEncodeToUInt(const char *s)
1895 if ((s[0] == 'x' || s[0] == 'X') && s[1] == 0)
1896 return kAtomTypeWildcard;
1897 if (s[0] >= '0' && s[0] <= '9')
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};
1906 if (*s == ',' || *s == '.' || *s == ';' || *s == ':') {
1907 /* Variant (only [0-9a-z] are allowed) */
1909 if (*s >= '0' && *s <= '9')
1911 else if (*s >= 'A' && *s <= 'Z')
1913 else if (*s >= 'a' && *s <= 'z')
1915 else n = (*s % 36) + 1; /* Map to something non-zero */
1916 t += n * (96*96*96*96);
1919 n = (*s - 0x20) % 96; /* Map to 1..95 */
1920 if (i == 0 && n < 32)
1922 t += n * s_coeff[i];
1928 AtomTypeDecodeToString(UInt type, char *s)
1930 static const UInt s_coeff[4] = {96*96*96, 96*96, 96, 1};
1938 if (type == kAtomTypeWildcard) {
1943 if (type < kAtomTypeMinimum) {
1944 snprintf(s, 6, "%d", type);
1947 for (i = 0; i < 4; i++) {
1948 s[i] = (type / s_coeff[i]) % 96;
1954 if ((variant = (type / (96*96*96*96)) % 96) != 0) {
1956 s[n + 1] = (variant <= 10 ? '0' + variant - 1 : 'a' + variant - 11);
1963 ElementToInt(const char *s)
1967 for (i = 0, p = gDispAtomParameters; i < gCountDispAtomParameters; i++, p++) {
1968 if (p->name[0] == toupper(s[0]) && p->name[1] == tolower(s[1]))
1975 ElementToString(Int elem, char *s)
1977 if (elem >= 0 && elem < gCountDispAtomParameters) {
1978 const char *cs = gDispAtomParameters[elem].name;
1987 AtomNameToElement(const char *s)
1991 /* $element = ($name =~ /([A-Za-z]{1,2})/); # in Perl */
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]);
2005 return ElementToInt(element);
2009 GuessAtomicNumber(const char *name, Double weight)
2015 for (i = 0, p = gDispAtomParameters; i < gCountDispAtomParameters; i++, p++) {
2016 if (p->weight > 0.0 && fabs(weight - p->weight) < 0.1)
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]);
2032 return ElementToInt(buf);
2036 WeightForAtomicNumber(Int elem)
2038 if (elem >= 1 && elem < gCountDispAtomParameters)
2039 return gDispAtomParameters[elem].weight;
2044 RadiusForAtomicNumber(Int elem)
2046 if (elem >= 1 && elem < gCountDispAtomParameters)
2047 return gDispAtomParameters[elem].radius;