MolLib/MainView.c

   1 /*
   2  *  MainView.c
   3  *
   4  *  Created by Toshi Nagata on 06/07/30.
   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 /*  On MinGW, #include of <GL/gl.h> should be before something in MolLib.h.
  18  Therefore, we include MainView.h first separately and then include
  19  the rest of MolLib.h.  */
  20 #include "MainView.h"
  21 #include "MolLib.h"
  22 #include "Ruby_bind/Molby_extern.h"
  23
  24 #include "MD/MDCore.h"
  25 #include "MD/MDGraphite.h"
  26
  27 #include <stdlib.h>
  28 #include <string.h>
  29 #include <ctype.h>
  30
  31 #define biso2radius(r) ((r) > 0.5 ? sqrt((r) / 78.9568352087147) : 0.08)
  32
  33 /*  pp is a pointer to Parameter record  */
  34 #define VDW_RADIUS(pp) (pp->vdw_radius == 0.0 ? pp->radius * 0.51 + 1.20 : pp->vdw_radius)
  35
  36 /*  Invalid bond/angle/torsion value, used in internal cache  */
  37 const Double kInvalidValue = -10000000.0;
  38
  39 #pragma mark ==== MainView public methods ====
  40
  41 void
  42 MainView_setViewObject(MainView *mview, void *ref)
  43 {
  44         static GLdouble sIdentity[16] = {
  45                 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1
  46         };
  47         if (mview != NULL) {
  48                 if (mview->ref == NULL) {
  49                         /*  Initialize GUI-only members  */
  50                         mview->mode = kTrackballRotateMode;
  51                         mview->tempAtoms[0] = mview->tempAtoms[1] = -1;
  52                         memmove(mview->modelview_matrix, sIdentity, sizeof(sIdentity));
  53                         memmove(mview->projection_matrix, sIdentity, sizeof(sIdentity));
  54                         mview->tableCache = IntGroupNew();
  55                         mview->tableSelection = IntGroupNew();
  56                 }
  57                 mview->ref = ref;  /*  No retain  */
  58                 IntGroupClear(mview->tableCache);
  59                 IntGroupClear(mview->tableSelection);
  60                 MoleculeCallback_notifyModification(mview->mol, 0);
  61         }
  62 }
  63
  64 int
  65 MainView_isAtomHidden(MainView *mview, int index)
  66 {
  67         if (index < 0 || index >= mview->mol->natoms)
  68                 return 1;
  69         else if (mview->visibleFlags == NULL)
  70                 return 0;
  71         else return (mview->visibleFlags[index] == 0);
  72 /*      Atom *ap = ATOM_AT_INDEX(mview->mol->atoms, index);
  73         if (ap->exflags & kAtomHiddenFlag)
  74                 return 1;
  75         if (!mview->showHydrogens && ap->atomicNumber == 1)
  76                 return 1;
  77         if (!mview->showDummyAtoms && ap->atomicNumber == 0)
  78                 return 1;
  79         return 0; */
  80 }
  81
  82 void
  83 MainView_refreshCachedInfo(MainView *mview)
  84 {
  85         Molecule *mol;
  86         int i, j, n1, n2, n3, n4;
  87         Byte f1, f2, f3, f4;
  88         Atom *ap;
  89
  90         if (mview == NULL || (mol = mview->mol) == NULL)
  91                 return;
  92
  93         /*  Rebuild internal caches  */
  94
  95         /*  Visible flags  */
  96         if (mview->visibleFlags == NULL)
  97                 mview->visibleFlags = (Byte *)malloc(mol->natoms);
  98         else
  99                 mview->visibleFlags = (Byte *)realloc(mview->visibleFlags, mol->natoms);
 100         memset(mview->visibleFlags, 0, mol->natoms);
 101         mview->countHidden = 0;
 102         for (i = 0, ap = mol->atoms; i < mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 103                 if (ap->exflags & kAtomHiddenFlag) {
 104                         mview->countHidden++;
 105                         continue;
 106                 }
 107                 if (!mview->showHydrogens && ap->atomicNumber == 1)
 108                         continue;
 109                 if (!mview->showDummyAtoms && ap->atomicNumber == 0)
 110                         continue;
 111                 mview->visibleFlags[i] = 1;
 112         }
 113
 114         /*  Selection (only temporary used within this function)  */
 115         {
 116                 IntGroup *sel = MoleculeGetSelection(mol);
 117                 if (sel != NULL && IntGroupGetCount(sel) > 0) {
 118                         for (i = 0; (n1 = IntGroupGetStartPoint(sel, i)) >= 0; i++) {
 119                                 n2 = IntGroupGetEndPoint(sel, i);
 120                                 if (n2 > mol->natoms)
 121                                         n2 = mol->natoms;
 122                                 while (n1 < n2)
 123                                         mview->visibleFlags[n1++] |= 2;
 124                         }
 125                 }
 126         }
 127
 128         /*  Table caches  */
 129         if (mview->tableCache == NULL)
 130                 mview->tableCache = IntGroupNew();
 131         if (mview->tableSelection == NULL)
 132                 mview->tableSelection = IntGroupNew();
 133         IntGroupClear(mview->tableCache);
 134         IntGroupClear(mview->tableSelection);
 135
 136         if (mview->tableIndex == kMainViewAtomTableIndex ||
 137                 mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Atoms */
 138                 for (i = j = 0; i < mol->natoms; i++) {
 139                         f1 = mview->visibleFlags[i];
 140                         if ((f1 & 1) != 0) {
 141                                 IntGroupAdd(mview->tableCache, i, 1);
 142                                 if ((f1 & 2) != 0)
 143                                         IntGroupAdd(mview->tableSelection, j, 1);
 144                                 j++;
 145                         }
 146                 }
 147         } else if (mview->tableIndex == kMainViewBondTableIndex) {  /* Bonds */
 148                 for (i = j = 0; i < mol->nbonds; i++) {
 149                         n1 = mol->bonds[i * 2];
 150                         n2 = mol->bonds[i * 2 + 1];
 151                         f1 = mview->visibleFlags[n1];
 152                         f2 = mview->visibleFlags[n2];
 153                         if ((f1 & 1) != 0 && (f2 & 1) != 0) {
 154                                 IntGroupAdd(mview->tableCache, i, 1);
 155                                 if ((f1 & 2) != 0 && (f2 & 2) != 0)
 156                                         IntGroupAdd(mview->tableSelection, j, 1);
 157                                 j++;
 158                         }
 159                 }
 160         } else if (mview->tableIndex == kMainViewAngleTableIndex) {  /* Angles */
 161                 for (i = j = 0; i < mol->nangles; i++) {
 162                         n1 = mol->angles[i * 3];
 163                         n2 = mol->angles[i * 3 + 1];
 164                         n3 = mol->angles[i * 3 + 2];
 165                         f1 = mview->visibleFlags[n1];
 166                         f2 = mview->visibleFlags[n2];
 167                         f3 = mview->visibleFlags[n3];
 168                         if ((f1 & 1) != 0 && (f2 & 1) != 0 && (f3 & 1) != 0) {
 169                                 IntGroupAdd(mview->tableCache, i, 1);
 170                                 if ((f1 & 2) != 0 && (f2 & 2) != 0 && (f3 & 2) != 0)
 171                                         IntGroupAdd(mview->tableSelection, j, 1);
 172                                 j++;
 173                         }
 174                 }
 175         } else if (mview->tableIndex == kMainViewDihedralTableIndex) {  /* Dihedrals */
 176                 for (i = j = 0; i < mol->ndihedrals; i++) {
 177                         n1 = mol->dihedrals[i * 4];
 178                         n2 = mol->dihedrals[i * 4 + 1];
 179                         n3 = mol->dihedrals[i * 4 + 2];
 180                         n4 = mol->dihedrals[i * 4 + 3];
 181                         f1 = mview->visibleFlags[n1];
 182                         f2 = mview->visibleFlags[n2];
 183                         f3 = mview->visibleFlags[n3];
 184                         f4 = mview->visibleFlags[n4];
 185                         if ((f1 & 1) != 0 && (f2 & 1) != 0 && (f3 & 1) != 0 && (f4 & 1) != 0) {
 186                                 IntGroupAdd(mview->tableCache, i, 1);
 187                                 if ((f1 & 2) != 0 && (f2 & 2) != 0 && (f3 & 2) != 0 && (f4 & 2) != 0)
 188                                         IntGroupAdd(mview->tableSelection, j, 1);
 189                                 j++;
 190                         }
 191                 }
 192         } else if (mview->tableIndex == kMainViewImproperTableIndex) {  /* Impropers */
 193                 for (i = j = 0; i < mol->nimpropers; i++) {
 194                         n1 = mol->impropers[i * 4];
 195                         n2 = mol->impropers[i * 4 + 1];
 196                         n3 = mol->impropers[i * 4 + 2];
 197                         n4 = mol->impropers[i * 4 + 3];
 198                         f1 = mview->visibleFlags[n1];
 199                         f2 = mview->visibleFlags[n2];
 200                         f3 = mview->visibleFlags[n3];
 201                         f4 = mview->visibleFlags[n4];
 202                         if ((f1 & 1) != 0 && (f2 & 1) != 0 && (f3 & 1) != 0 && (f4 & 1) != 0) {
 203                                 IntGroupAdd(mview->tableCache, i, 1);
 204                                 if ((f1 & 2) != 0 && (f2 & 2) != 0 && (f3 & 2) != 0 && (f4 & 2) != 0)
 205                                         IntGroupAdd(mview->tableSelection, j, 1);
 206                                 j++;
 207                         }
 208                 }
 209         } else {
 210                 /*  Cache is left empty (not used)  */
 211         }
 212
 213 //      } else if (mview->tableIndex == kMainViewParameterTableIndex ||  /* Parameter infos */
 214 //                         mview->tableIndex == kMainViewUnitCellTableIndex) {   /* Unit cell infos */
 215 //              /*  Do nothing (tableCache will not be used)  */
 216 //      } else if (mview->tableIndex == kMainViewMOTableIndex) {  /* MO infos  */
 217 //              /*  Really no need to cache info, but create it anyway to simplify code  */
 218 //              if (mol->bset != NULL && mol->bset->ncomps > 0)
 219 //                      IntGroupAdd(mview->tableCache, 0, mol->bset->ncomps);
 220 //      }
 221
 222         /*  Clear internal selection flag  */
 223         for (i = 0; i < mol->natoms; i++) {
 224                 mview->visibleFlags[i] &= ~2;
 225         }
 226 }
 227
 228 #pragma mark ====== 2D/3D transform operations ======
 229
 230 void
 231 MainView_resizeToFit(MainView *mview)
 232 {
 233         Vector p;
 234         double f[4];
 235         if (mview == NULL || mview->mol == NULL)
 236                 return;
 237         if (mview->mol->natoms == 0) {
 238                 TrackballReset(mview->track);
 239                 return;
 240         }
 241         MoleculeCenterOfMass(mview->mol, &p, NULL);
 242         f[0] = -p.x / mview->dimension;
 243         f[1] = -p.y / mview->dimension;
 244         f[2] = -p.z / mview->dimension;
 245         TrackballSetTranslate(mview->track, f);
 246
 247         /*  Set scale
 248                 r0: the longest distance from the center of mass
 249                 r0 > dimension: scale = -log(r0/dimension)/log(10)  (negative)
 250                 r0 < dimension: scale = -log(atan2(r0, dimension*cot(15deg))*180deg/pi*2/30deg)/log(10) (positive)
 251         */
 252         {
 253                 int i;
 254                 Vector q;
 255                 Atom *ap;
 256                 double r0 = 0.0, r1, scale;
 257                 for (i = 0, ap = mview->mol->atoms; i < mview->mol->natoms; i++, ap = ATOM_NEXT(ap)) {
 258                         q = ap->r;
 259                         VecDec(q, p);
 260                         r1 = VecLength(q);
 261                         if (r1 > r0)
 262                                 r0 = r1;
 263                 }
 264                 r0 /= mview->dimension;
 265                 if (r0 < 1e-6)
 266                         scale = 0.0;
 267                 else if (r0 < 1.0)
 268                         scale = -log(atan2(r0, kCot15Deg) * kRad2Deg * 2 / 30.0) / kLog10;
 269                 else
 270                         scale = -log(r0) / kLog10;
 271                 TrackballSetScale(mview->track, scale);
 272         }
 273
 274         MainViewCallback_setNeedsDisplay(mview, 1);
 275 }
 276
 277 int
 278 MainView_convertScreenPositionToObjectPosition(MainView *mview, const GLfloat *screenPos, double *objectPos)
 279 {
 280         float rect[4];
 281     GLint viewport[4], n;
 282     GLfloat winZ;
 283     GLdouble posX, posY, posZ;
 284         if (mview == NULL)
 285                 return 0;
 286         MainViewCallback_frame(mview, rect);
 287     viewport[0] = viewport[1] = 0;
 288     viewport[2] = (GLint)(rect[2] - rect[0]);
 289     viewport[3] = (GLint)(rect[3] - rect[1]);
 290         MainViewCallback_lockFocus(mview);
 291     if (screenPos[2] >= 0.0)
 292         winZ = screenPos[2];
 293     else
 294         glReadPixels(screenPos[0], screenPos[1], 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
 295     gluUnProject(screenPos[0], screenPos[1], winZ, mview->modelview_matrix, mview->projection_matrix, viewport, &posX, &posY, &posZ);
 296     n = glGetError();
 297         MainViewCallback_unlockFocus(mview);
 298         objectPos[0] = posX;
 299         objectPos[1] = posY;
 300         objectPos[2] = posZ;
 301     if (n != GL_NO_ERROR || winZ == 1.0)
 302         return 0;
 303     else
 304         return 1;
 305 }
 306
 307 int
 308 MainView_convertObjectPositionToScreenPosition(MainView *mview, const double *objectPos, GLfloat *screenPos)
 309 {
 310         float rect[4];
 311     GLint viewport[4];
 312     GLdouble objX, objY, objZ;
 313         if (mview == NULL)
 314                 return 0;
 315         MainViewCallback_frame(mview, rect);
 316     viewport[0] = viewport[1] = 0;
 317     viewport[2] = (GLint)(rect[2] - rect[0]);
 318     viewport[3] = (GLint)(rect[3] - rect[1]);
 319     gluProject(objectPos[0], objectPos[1], objectPos[2], mview->modelview_matrix, mview->projection_matrix, viewport, &objX, &objY, &objZ);
 320     if (glGetError() == GL_NO_ERROR) {
 321                 screenPos[0] = objX;
 322                 screenPos[1] = objY;
 323                 screenPos[2] = objZ;
 324         /*      fprintf(stderr, "object(%.3f,%.3f,%.3f) screen(%.3f,%.3f,%.3f)\n", objectPos[0], objectPos[1], objectPos[2], screenPos[0], screenPos[1], screenPos[2]); */
 325                 return 1;
 326         } else return 0;
 327 }
 328
 329 int
 330 MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex1, int *outIndex2, int mouseCheck, int ignoreExpandedAtoms)
 331 {
 332         float screenPos[3];
 333         double op[3], oq[3], pqlen, pqlen2;
 334         Vector pq, pa, v1, r1, r2;
 335     int i, natoms, nbonds;
 336         float r, d2, z;
 337         const Atom *ap, *bp;
 338         const ElementPar *ep;
 339         const int *ip;
 340         Molecule *mol;
 341         float minDepth;
 342         int n1, n2;
 343
 344         if (mview == NULL || mview->mol == NULL) {
 345                 *outIndex1 = *outIndex2 = -1;
 346                 return 0;
 347         }
 348         mol = mview->mol;
 349
 350         screenPos[0] = mousePos[0];
 351         screenPos[1] = mousePos[1];
 352         screenPos[2] = -1.0;
 353         if (MainView_convertScreenPositionToObjectPosition(mview, screenPos, op) == 0)
 354                 return 0;  /*  Nothing is here  */
 355
 356         /*  PQ is the part of the eyesight line in the visible area  */
 357         screenPos[2] = 0.0;
 358         MainView_convertScreenPositionToObjectPosition(mview, screenPos, op);
 359         screenPos[2] = 1.0;
 360         MainView_convertScreenPositionToObjectPosition(mview, screenPos, oq);
 361         pq.x = oq[0] - op[0];
 362         pq.y = oq[1] - op[1];
 363         pq.z = oq[2] - op[2];
 364         pqlen2 = VecLength2(pq);
 365         pqlen = sqrt(pqlen2);
 366         natoms = mol->natoms;
 367         n1 = n2 = -1;
 368         minDepth = 100.0;
 369         for (i = 0; i < natoms; i++) {
 370                 Vector pq1, pa1;
 371                 double pq1len2, pq1len;
 372                 if (mouseCheck && i % 50 == 0 && MainViewCallback_mouseCheck(mview))
 373                         return 0;  /*  If mouse event is detected return immediately  */
 374                 /*  Examine if an atom is visible or not  */
 375                 /*  The distance of the atom center (A) from line PQ: */
 376                 /*    d = |VecCross(PA,PQ)|/|PQ|  */
 377                 /*  z = VecDot(PA,PQ)/|PQ|^2 - sqrt(r^2 - d^2)/|PQ|  */
 378         ap = ATOM_AT_INDEX(mol->atoms, i);
 379                 if (ap == NULL)
 380                         continue;
 381                 if (MainView_isAtomHidden(mview, i))
 382                         continue;
 383                 if (ignoreExpandedAtoms && SYMOP_ALIVE(ap->symop))
 384                         continue;
 385                 r1 = ap->r;
 386         /*      if (mol->is_xtal_coord)
 387                         TransformVec(&r1, mol->cell->tr, &r1); */
 388                 pa.x = r1.x - op[0];
 389                 pa.y = r1.y - op[1];
 390                 pa.z = r1.z - op[2];
 391                 if (mview->showEllipsoids && ap->aniso != NULL) {
 392                         /*  Convert to ellipsoid principal axes  */
 393                         Mat33 m1;
 394                         Aniso *anp = ap->aniso;
 395                         MatrixInvert(m1, anp->pmat);
 396                         MatrixVec(&pq1, m1, &pq);
 397                         MatrixVec(&pa1, m1, &pa);
 398                         r = mview->probabilityScale;
 399                         pq1len2 = VecLength2(pq1);
 400                         pq1len = sqrt(pq1len2);
 401                 } else {
 402                         if (mview->showEllipsoids) {
 403                                 r = biso2radius(ap->tempFactor) * mview->probabilityScale;
 404                         } else {
 405                                 ep = &(gElementParameters[ap->atomicNumber]);
 406                                 if (ep == NULL)
 407                                         continue;
 408                                 r = VDW_RADIUS(ep) * mview->atomRadius;
 409                         }
 410                         pa1 = pa;
 411                         pq1 = pq;
 412                         pq1len2 = pqlen2;
 413                         pq1len = pqlen;
 414                 }
 415                 VecCross(v1, pa1, pq1);
 416                 d2 = VecLength2(v1) / pq1len2;
 417                 if (d2 > r * r)
 418                         continue;  /*  Not visible  */
 419                 z = VecDot(pa1, pq1) / pq1len2 - sqrt(r * r - d2) / pq1len;
 420                 if (z < 0.0 || z > 1.0)
 421                         continue;  /*  Out of viewing volume  */
 422                 if (z < minDepth) {
 423                         minDepth = z;
 424                         n1 = i;
 425                 }
 426         }
 427         nbonds = mol->nbonds;
 428         for (i = 0; i < nbonds; i++) {
 429                 Vector vx, vy, vz, vv, vp;
 430                 Double wb, wa, t, wx, blen;
 431                 if (mouseCheck && i % 50 == 0 && MainViewCallback_mouseCheck(mview))
 432                         return 0;  /*  If mouse event is detected return immediately  */
 433                 /*  Examine if a bond is visible or not  */
 434                 ip = &(mol->bonds[i * 2]);
 435                 ap = ATOM_AT_INDEX(mol->atoms, ip[0]);
 436                 bp = ATOM_AT_INDEX(mol->atoms, ip[1]);
 437                 if (MainView_isAtomHidden(mview, ip[0]) || MainView_isAtomHidden(mview, ip[1]))
 438                         continue;
 439                 if (ignoreExpandedAtoms && SYMOP_ALIVE(ap->symop) && SYMOP_ALIVE(bp->symop))
 440                         continue;
 441                 /*  vx/vy/vz is a orthonormal base in which AB parallels the x-axis
 442                     and AP in in the xy plane  */
 443                 /*  vp and vv is AP and PQ in that coordinate system  */
 444                 r1 = ap->r;
 445                 r2 = bp->r;
 446         /*      if (mol->is_xtal_coord) {
 447                         TransformVec(&r1, mol->cell->tr, &r1);
 448                         TransformVec(&r2, mol->cell->tr, &r2);
 449                 } */
 450                 v1.x = op[0] - r1.x;
 451                 v1.y = op[1] - r1.y;
 452                 v1.z = op[2] - r1.z;
 453                 VecSub(vx, r2, r1);
 454                 blen = sqrt(VecLength2(vx));
 455                 if (blen < 1e-10)
 456                         continue;
 457                 vx.x /= blen;
 458                 vx.y /= blen;
 459                 vx.z /= blen;
 460                 VecCross(vz, vx, v1);
 461                 if (NormalizeVec(&vz, &vz))
 462                         continue;
 463                 VecCross(vy, vz, vx);
 464                 vp.x = VecDot(v1, vx);
 465                 vp.y = VecDot(v1, vy);
 466                 vp.z = VecDot(v1, vz);
 467                 vv.x = VecDot(pq, vx);
 468                 vv.y = VecDot(pq, vy);
 469                 vv.z = VecDot(pq, vz);
 470                 /*  The bond surface is y^2 + z^2 = r^2, 0 <= x <= 1  */
 471                 /*  The eyesight line is (x,y,z) = vp + t * vv, 0 <= t <= 1  */
 472                 /*  The crossing point: t = (-vv.y*vp.y - sqrt((vv.y*vp.y)^2 - (vv.y^2+vv.z^2)(vp.y^2-r^2)))/(vv.y^2+vv.z^2)  */
 473                 /*  (Note that vp.z = 0 by definition)  */
 474                 r = mview->bondRadius;
 475                 wb = vv.y * vp.y;
 476                 wa = vv.y * vv.y + vv.z * vv.z;
 477                 d2 = wb * wb - wa * (vp.y * vp.y - r * r);
 478                 if (d2 < 0)
 479                         continue;  /*  Not visible  */
 480                 t = (-wb - sqrt(d2)) / wa;
 481                 if (t < 0 || t > 1)
 482                         continue;  /*  Out of visible volume  */
 483                 wx = vp.x + t * vv.x;
 484                 if (wx < 0 || wx > blen)
 485                         continue;  /*  Outside of bond  */
 486                 if (t < minDepth) {
 487                         minDepth = t;
 488                         n1 = ip[0];
 489                         n2 = ip[1];
 490                 }
 491         }
 492         *outIndex1 = n1;
 493         *outIndex2 = n2;
 494         return (n1 >= 0 || n2 >= 0);
 495 }
 496
 497 int
 498 MainView_screenCenterPointOfAtom(MainView *mview, int index, float *outScreenPos)
 499 {
 500         const Atom *ap;
 501         const ElementPar *dp;
 502         Vector cv, pv, v;
 503         Double rad, w;
 504         double p[3];
 505
 506         if (mview == NULL || mview->mol == NULL || index < 0 || index >= mview->mol->natoms)
 507                 return 0;
 508         ap = ATOM_AT_INDEX(mview->mol->atoms, index);
 509
 510         /*  Camera position in object coordinates  */
 511         cv = mview->camera;
 512
 513         /*  The atom position (in Cartesian)  */
 514         v = ap->r;
 515 /*      if (mview->mol->is_xtal_coord)
 516                 TransformVec(&v, mview->mol->cell->tr, &v); */
 517
 518         /*  The vector from atom center to camera  */
 519         VecSub(pv, cv, v);
 520
 521         /*  Get the surface point of the ellipsoid/sphere along the camera vector  */
 522         if (mview->showEllipsoids) {
 523                 Mat33 m1;
 524                 Aniso *anp = ap->aniso;
 525                 if (anp != NULL) {
 526                         Vector vx, vy, vz;
 527                         MatrixInvert(m1, anp->pmat);
 528                         /*  Convert the 'screen plane' vectors to ellipsoid principal axes  */
 529                         vy = mview->up;
 530                         VecCross(vx, mview->lookto, vy);
 531                         MatrixVec(&vx, m1, &vx);
 532                         MatrixVec(&vy, m1, &vy);
 533                         /*  Touching point of the 'screen plane' to the ellipsoid  */
 534                         VecCross(vz, vx, vy);
 535                         w = mview->probabilityScale / VecLength(vz) * 1.1;
 536                         VecScaleSelf(vz, w);
 537                         MatrixVec(&vz, anp->pmat, &vz);
 538                         /*  The crossing point of the camera vector with the touching plane */
 539                         w = fabs(VecDot(pv, vz) / VecLength2(pv));
 540                         VecScaleSelf(pv, w);
 541                 } else {
 542                         w = mview->probabilityScale * biso2radius(ap->tempFactor) / VecLength(pv) * 1.1;
 543                         VecScaleSelf(pv, w);
 544                 }
 545         } else {
 546                 dp = &(gElementParameters[ap->atomicNumber]);
 547                 rad = VDW_RADIUS(dp) * mview->atomRadius;
 548                 w = rad / VecLength(pv) * 1.1;
 549                 VecScaleSelf(pv, w);
 550         }
 551         VecInc(v, pv);
 552
 553         /*  Transform to screen coordinate  */
 554         p[0] = v.x;
 555         p[1] = v.y;
 556         p[2] = v.z;
 557         return MainView_convertObjectPositionToScreenPosition(mview, p, outScreenPos);
 558 }
 559
 560 void
 561 MainView_getCamera(MainView *mview, Vector *outCamera, Vector *outLookAt, Vector *outUp)
 562 {
 563         if (mview != NULL) {
 564                 *outCamera = mview->camera;
 565                 *outLookAt = mview->lookat;
 566                 *outUp = mview->up;
 567         }
 568 }
 569
 570 #pragma mark ====== Draw model ======
 571
 572 static void
 573 enableLighting(void)
 574 {
 575         static GLfloat pos[] = {0.0f, 0.0f, 1.0f, 0.0f};
 576         glEnable(GL_LIGHTING);
 577 #if __WXMAC__
 578         /*  On Mac OS 10.6, redefining the light position seems to be necessary
 579                 every time lighting is made enabled  */
 580         glMatrixMode(GL_MODELVIEW);
 581         glPushMatrix();
 582         glLoadIdentity();
 583         glLightfv(GL_LIGHT0, GL_POSITION, pos);
 584         glPopMatrix();
 585 #endif
 586 }
 587
 588 void
 589 MainView_initializeOpenGL(void)
 590 {
 591         static GLfloat ambient[] = {0.6, 0.6, 0.6, 1.0};  // Some white ambient light.
 592         static GLfloat diffuse[] = {1.0, 1.0, 1.0, 1.0};  // A white light.
 593
 594         glEnable(GL_DEPTH_TEST);
 595         glEnable(GL_LIGHTING);
 596
 597         //  Set the ambient light
 598         glLightModelfv(GL_LIGHT_MODEL_AMBIENT, ambient);
 599
 600         //  Set the light and switch it on
 601         glLightfv(GL_LIGHT0, GL_DIFFUSE, diffuse);
 602         glEnable(GL_LIGHT0);
 603
 604         //  Enable blending
 605         glEnable(GL_BLEND);
 606         glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
 607 }
 608
 609 /*  Get orthogonal unit vectors  */
 610 static int
 611 getOrthogonalVectors(const GLfloat *ap, GLfloat *bp, GLfloat *cp)
 612 {
 613     Double ra, rb;
 614     ra = sqrt(ap[0] * ap[0] + ap[1] * ap[1] + ap[2] * ap[2]);
 615     if (ra < 1e-20)
 616         return 0;
 617     ra = 1 / ra;
 618     if (fabs(ap[0]) < fabs(ap[1])) {
 619         if (fabs(ap[0]) < fabs(ap[2])) {
 620             bp[0] = 0;
 621             bp[1] = -ap[2];
 622             bp[2] = ap[1];
 623         } else {
 624             bp[0] = ap[1];
 625             bp[1] = -ap[0];
 626             bp[2] = 0;
 627         }
 628     } else {
 629         if (fabs(ap[1]) < fabs(ap[2])) {
 630             bp[0] = -ap[2];
 631             bp[1] = 0;
 632             bp[2] = ap[0];
 633         } else {
 634             bp[0] = ap[1];
 635             bp[1] = -ap[0];
 636             bp[2] = 0;
 637         }
 638     }
 639     rb = 1 / sqrt(bp[0] * bp[0] + bp[1] * bp[1] + bp[2] * bp[2]);
 640     bp[0] *= rb;
 641     bp[1] *= rb;
 642     bp[2] *= rb;
 643     cp[0] = ra * (ap[1] * bp[2] - ap[2] * bp[1]);
 644     cp[1] = ra * (ap[2] * bp[0] - ap[0] * bp[2]);
 645     cp[2] = ra * (ap[0] * bp[1] - ap[1] * bp[0]);
 646 /*    printf("a = (%f, %f, %f) b = (%f, %f, %f) c = (%f, %f, %f)\n",
 647         ap[0], ap[1], ap[2], bp[0], bp[1], bp[2], cp[0], cp[1], cp[2]); */
 648     return 1;
 649 }
 650
 651 static GLfloat sSinCache[81];
 652 static int sSinCacheSect = 0;
 653
 654 static int
 655 setSinCache(int sect)
 656 {
 657     int n, m, i;
 658     m = sect / 4;
 659     n = m * 4;
 660     if (n >= 64)
 661         n = 64;
 662     if (n != sSinCacheSect) {
 663         sSinCacheSect = n;
 664         for (i = 0; i <= m * 5; i++)
 665             sSinCache[i] = sin(3.14159265358979 * 2 / n * i);
 666     }
 667     return n;
 668 }
 669
 670 static void
 671 drawCylinder(const GLfloat *a, const GLfloat *b, GLfloat r, int sect, int closed)
 672 {
 673     GLfloat *c, *s;
 674     int n, i;
 675         float nx, ny, nz;
 676     GLfloat d[3], v[3], w[3];
 677     n = setSinCache(sect);
 678     if (n <= 0)
 679         return;
 680     s = sSinCache;
 681     c = &sSinCache[n/4];
 682     d[0] = b[0] - a[0];
 683     d[1] = b[1] - a[1];
 684     d[2] = b[2] - a[2];
 685     if (getOrthogonalVectors(d, v, w) == 0)
 686         return;
 687     glBegin(GL_QUAD_STRIP);
 688     for (i = 0; i <= n; i++) {
 689         nx = v[0] * c[i] + w[0] * s[i];
 690         ny = v[1] * c[i] + w[1] * s[i];
 691         nz = v[2] * c[i] + w[2] * s[i];
 692         glNormal3f(nx, ny, nz);
 693         glVertex3f(a[0] + r * nx, a[1] + r * ny, a[2] + r * nz);
 694         glVertex3f(b[0] + r * nx, b[1] + r * ny, b[2] + r * nz);
 695     }
 696     glEnd();
 697         if (closed) {
 698                 glBegin(GL_TRIANGLE_FAN);
 699                 glNormal3f(-d[0], -d[1], -d[2]);
 700                 for (i = 0; i <= n; i++) {
 701                         nx = v[0] * c[i] + w[0] * s[i];
 702                         ny = v[1] * c[i] + w[1] * s[i];
 703                         nz = v[2] * c[i] + w[2] * s[i];
 704                         glVertex3f(a[0] + r * nx, a[1] + r * ny, a[2] + r * nz);
 705                 }
 706                 glEnd();
 707                 glBegin(GL_TRIANGLE_FAN);
 708                 glNormal3f(d[0], d[1], d[2]);
 709                 for (i = 0; i <= n; i++) {
 710                         nx = v[0] * c[i] + w[0] * s[i];
 711                         ny = v[1] * c[i] + w[1] * s[i];
 712                         nz = v[2] * c[i] + w[2] * s[i];
 713                         glVertex3f(b[0] + r * nx, b[1] + r * ny, b[2] + r * nz);
 714                 }
 715                 glEnd();
 716         }
 717 }
 718
 719 static void
 720 drawCone(const GLfloat *a, const GLfloat *b, GLfloat r, int sect, int closed)
 721 {
 722     GLfloat *c, *s;
 723     int n, i;
 724     GLfloat d[3], v[3], w[3];
 725         Vector p1, nv;
 726     n = setSinCache(sect);
 727     if (n <= 0)
 728         return;
 729     s = sSinCache;
 730     c = &sSinCache[n/4];
 731     d[0] = b[0] - a[0];
 732     d[1] = b[1] - a[1];
 733     d[2] = b[2] - a[2];
 734     if (getOrthogonalVectors(d, v, w) == 0)
 735         return;
 736     glBegin(GL_TRIANGLE_FAN);
 737         nv.x = d[0];
 738         nv.y = d[1];
 739         nv.z = d[2];
 740         NormalizeVec(&nv, &nv);
 741         glNormal3f(nv.x, nv.y, nv.z);
 742         glVertex3f(b[0], b[1], b[2]);
 743     for (i = 0; i <= n; i++) {
 744         nv.x = v[0] * c[i] + w[0] * s[i];
 745         nv.y = v[1] * c[i] + w[1] * s[i];
 746         nv.z = v[2] * c[i] + w[2] * s[i];
 747                 glNormal3f(nv.x, nv.y, nv.z);
 748                 p1.x = a[0] + r * nv.x;
 749                 p1.y = a[1] + r * nv.y;
 750                 p1.z = a[2] + r * nv.z;
 751         glNormal3f(nv.x, nv.y, nv.z);
 752                 glVertex3f(p1.x, p1.y, p1.z);
 753     }
 754     glEnd();
 755         if (closed) {
 756                 glBegin(GL_TRIANGLE_FAN);
 757                 glNormal3f(d[0], d[1], d[2]);
 758                 for (i = 0; i <= n; i++) {
 759                         p1.x = a[0] + r * (v[0] * c[i] + w[0] * s[i]);
 760                         p1.y = a[1] + r * (v[1] * c[i] + w[1] * s[i]);
 761                         p1.z = a[2] + r * (v[2] * c[i] + w[2] * s[i]);
 762                         glVertex3f(p1.x, p1.y, p1.z);
 763                 }
 764                 glEnd();
 765         }
 766 }
 767
 768 static void
 769 drawSphere(const GLfloat *p, GLfloat r, int sect)
 770 {
 771     GLfloat *c, *s;
 772     int n, i, j;
 773     n = setSinCache(sect);
 774     if (n <= 0)
 775         return;
 776     s = sSinCache;
 777     c = &sSinCache[n/4];
 778     for (i = 0; i <= n; i++) {
 779         glBegin(GL_QUAD_STRIP);
 780         for (j = 1; j <= n / 2 - 1; j++) {
 781             glNormal3f(s[j] * c[i], s[j] * s[i], c[j]);
 782             glVertex3f(r * s[j] * c[i] + p[0], r * s[j] * s[i] + p[1], r * c[j] + p[2]);
 783             glNormal3f(s[j] * c[i+1], s[j] * s[i+1], c[j]);
 784             glVertex3f(r * s[j] * c[i+1] + p[0], r * s[j] * s[i+1] + p[1], r * c[j] + p[2]);
 785         }
 786         glEnd();
 787     }
 788     glBegin(GL_TRIANGLE_FAN);
 789     glNormal3f(0, 0, 1);
 790     glVertex3f(p[0], p[1], r + p[2]);
 791     for (i = n; i >= 0; i--) {
 792         glNormal3f(s[1] * c[i], s[1] * s[i], c[1]);
 793         glVertex3f(r * s[1] * c[i] + p[0], r * s[1] * s[i] + p[1], r * c[1] + p[2]);
 794     }
 795     glEnd();
 796     glBegin(GL_TRIANGLE_FAN);
 797     glNormal3f(0, 0, -1);
 798     glVertex3f(p[0], p[1], -r + p[2]);
 799     for (i = 0; i <= n; i++) {
 800         glNormal3f(s[1] * c[i], s[1] * s[i], -c[1]);
 801         glVertex3f(r * s[1] * c[i] + p[0], r * s[1] * s[i] + p[1], -r * c[1] + p[2]);
 802     }
 803     glEnd();
 804 }
 805
 806 static void
 807 drawEllipsoid(const GLfloat *p, const GLfloat *v1, const GLfloat *v2, const GLfloat *v3, int sect)
 808 {
 809         GLfloat mat[16];
 810         static const GLfloat origin[3] = {0, 0, 0};
 811         mat[0] = v1[0]; mat[1] = v1[1]; mat[2] = v1[2]; mat[3] = 0;
 812         mat[4] = v2[0]; mat[5] = v2[1]; mat[6] = v2[2]; mat[7] = 0;
 813         mat[8] = v3[0]; mat[9] = v3[1]; mat[10] = v3[2]; mat[11] = 0;
 814         mat[12] = p[0]; mat[13] = p[1]; mat[14] = p[2]; mat[15] = 1;
 815     glMatrixMode(GL_MODELVIEW);
 816         glPushMatrix();
 817         glMultMatrixf(mat);
 818         glEnable(GL_NORMALIZE);
 819         //      glutSolidSphere(1, sect, sect); /* Is this faster than my code? */
 820         drawSphere(origin, 1, sect);
 821         glDisable(GL_NORMALIZE);
 822         glPopMatrix();
 823 }
 824
 825 static char *
 826 temporarySelection(MainView *mview, int flags, int clickCount, int ignoreExpandedAtoms)
 827 {
 828         char *selectFlags;
 829         int i, natoms;
 830         const Atom *ap;
 831         Double rect[4];
 832         natoms = mview->mol->natoms;
 833         selectFlags = (char *)calloc(sizeof(char), natoms);
 834         if (selectFlags == NULL)
 835                 return NULL;
 836         if (clickCount > 0) {
 837                 int n1, n2;
 838                 if (MainView_findObjectAtPoint(mview, mview->dragStartPos, &n1, &n2, 0, 0)) {
 839                         if (n1 >= 0 && n1 < natoms)
 840                                 selectFlags[n1] = 1;
 841                         if (n2 >= 0 && n2 < natoms)
 842                                 selectFlags[n2] = 1;
 843                 }
 844         } else {
 845                 if (mview->dragStartPos[0] < mview->dragEndPos[0]) {
 846                         rect[0] = mview->dragStartPos[0];
 847                         rect[2] = mview->dragEndPos[0];
 848                 } else {
 849                         rect[0] = mview->dragEndPos[0];
 850                         rect[2] = mview->dragStartPos[0];
 851                 }
 852                 if (mview->dragStartPos[1] < mview->dragEndPos[1]) {
 853                         rect[1] = mview->dragStartPos[1];
 854                         rect[3] = mview->dragEndPos[1];
 855                 } else {
 856                         rect[1] = mview->dragEndPos[1];
 857                         rect[3] = mview->dragStartPos[1];
 858                 }
 859                 for (i = 0; i < natoms; i++) {
 860                         ap = ATOM_AT_INDEX(mview->mol->atoms, i);
 861                         if (ap == NULL)
 862                                 continue;
 863                         if (MainView_isAtomHidden(mview, i))
 864                                 continue;
 865                         if (ignoreExpandedAtoms && SYMOP_ALIVE(ap->symop))
 866                                 continue;
 867                         if (mview->draggingMode == kMainViewSelectingRegion) {
 868                                 /*  Check if this atom is within the selection rectangle  */
 869                                 double objectPos[3];
 870                                 GLfloat screenPos[3];
 871                                 Vector r1;
 872                                 r1 = ap->r;
 873                         /*      if (mview->mol->is_xtal_coord)
 874                                         TransformVec(&r1, mview->mol->cell->tr, &r1); */
 875                                 objectPos[0] = r1.x;
 876                                 objectPos[1] = r1.y;
 877                                 objectPos[2] = r1.z;
 878                                 if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, screenPos) && screenPos[0] >= rect[0] && screenPos[0] <= rect[2] && screenPos[1] >= rect[1] && screenPos[1] <= rect[3])
 879                                         selectFlags[i] = 1;
 880                                 else selectFlags[i] = 0;
 881                         }
 882                 }
 883         }
 884         if (flags & kShiftKeyMask) {
 885                 for (i = 0; i < natoms; i++)
 886                         selectFlags[i] ^= (MoleculeIsAtomSelected(mview->mol, i) != 0);
 887         }
 888         return selectFlags;
 889 }
 890
 891 static void
 892 drawGraphite(MainView *mview)
 893 {
 894         static GLfloat sDarkCyanColor[] = {0, 0.75, 0.75, 1};
 895         MDArena *arena;
 896         MDGraphiteArena *graphite;
 897         Vector xaxis, yaxis, zaxis, origin;
 898         Double R;
 899         int i, j, i0, i1, j0, j1, ir;
 900         Double x, dx, y, dy, xx, yy;
 901         GLfloat p[12];
 902         if ((arena = mview->mol->arena) != NULL && (graphite = arena->graphite) != NULL) {
 903                 graphite_get_axes(graphite, &origin, &xaxis, &yaxis, &zaxis, &R);
 904         } else {
 905                 origin.x = origin.y = origin.z = 0.0;
 906                 xaxis.x = yaxis.y = zaxis.z = 1.0;
 907                 xaxis.y = xaxis.z = yaxis.x = yaxis.z = zaxis.x = zaxis.y = 0.0;
 908                 R = 1.42;
 909         }
 910         i0 = -(mview->showGraphite / 2) - 1;
 911         i1 = i0 + mview->showGraphite + 1;
 912         j0 = -(mview->showGraphite / 2);
 913         j1 = j0 + mview->showGraphite;
 914         dx = 0.5 * R;
 915         dy = 0.866025403784439 * R;
 916         glDisable(GL_LIGHTING);
 917         glColor3fv(sDarkCyanColor);
 918         for (i = i0; i <= i1; i++) {
 919                 for (j = j0; j <= j1; j++) {
 920                         Byte f1, f2, f3;
 921                         ir = (i % 2 == 0 ? 0 : 1);
 922                         x = 3 * i * dx;
 923                         y = (2 * j + ir) * dy;
 924                         yy = y - dy;
 925                         xx = x - 2 * dx;
 926                         p[0] = xaxis.x * xx + yaxis.x * y + origin.x;
 927                         p[1] = xaxis.y * xx + yaxis.y * y + origin.y;
 928                         p[2] = xaxis.z * xx + yaxis.z * y + origin.z;
 929                         xx += dx;
 930                         p[3] = xaxis.x * xx + yaxis.x * yy + origin.x;
 931                         p[4] = xaxis.y * xx + yaxis.y * yy + origin.y;
 932                         p[5] = xaxis.z * xx + yaxis.z * yy + origin.z;
 933                         xx += 2 * dx;
 934                         p[6] = xaxis.x * xx + yaxis.x * yy + origin.x;
 935                         p[7] = xaxis.y * xx + yaxis.y * yy + origin.y;
 936                         p[8] = xaxis.z * xx + yaxis.z * yy + origin.z;
 937                         xx += dx;
 938                         p[9] = xaxis.x * xx + yaxis.x * y + origin.x;
 939                         p[10] = xaxis.y * xx + yaxis.y * y + origin.y;
 940                         p[11] = xaxis.z * xx + yaxis.z * y + origin.z;
 941                         f1 = f2 = f3 = 1;
 942                         if (i == i0) {
 943                                 f1 = f2 = 0;
 944                                 if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
 945                                         continue;
 946                         } else if (i == i1) {
 947                                 f2 = f3 = 0;
 948                                 if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
 949                                         continue;
 950                         } else if (j == j1) {
 951                                 if (ir == 1) {
 952                                         f1 = f3 = 0;
 953                                 } else if (i == i0 + 1) {
 954                                         f1 = 0;
 955                                 } else if (i == i1 - 1) {
 956                                         f3 = 0;
 957                                 }
 958                         }
 959                         glBegin(GL_LINES);
 960                         if (f1) {
 961                                 glVertex3fv(p);
 962                                 glVertex3fv(p + 3);
 963                         }
 964                         if (f2) {
 965                                 glVertex3fv(p + 3);
 966                                 glVertex3fv(p + 6);
 967                         }
 968                         if (f3) {
 969                                 glVertex3fv(p + 6);
 970                                 glVertex3fv(p + 9);
 971                         }
 972                         glEnd();
 973                 }
 974         }
 975         enableLighting();
 976 }
 977
 978 static GLfloat sRedColor[] = {1, 0, 0, 1};
 979
 980 static void
 981 drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const Vector *periodicOffset)
 982 {
 983         const Atom *ap;
 984         const ElementPar *dp;
 985         int an1;
 986         int expanded = 0;
 987         Vector r1;
 988         GLfloat p[6];
 989         double pp[3];
 990         double rad;
 991         char label[16];
 992         GLfloat rgba[4];
 993         Transform *trp = NULL;
 994         int natoms = mview->mol->natoms;
 995         if (i1 >= natoms) {
 996                 /*  Extra 2 atoms for the bond being newly created  */
 997                 if (mview->draggingMode != kMainViewCreatingBond)
 998                         return;
 999                 if (mview->tempAtoms[i1 - natoms] >= 0)
1000                         return;  /*  Already drawn  */
1001                 ap = NULL;
1002                 an1 = 6;
1003                 r1 = mview->tempAtomPos[i1 - natoms];
1004                 label[0] = 0;
1005         } else {
1006                 ap = ATOM_AT_INDEX(mview->mol->atoms, i1);
1007                 if (ap == NULL)
1008                         return;
1009                 an1 = ap->atomicNumber;
1010                 r1 = ap->r;
1011                 strncpy(label, ap->aname, 4);
1012                 label[4] = 0;
1013                 if (SYMOP_ALIVE(ap->symop))
1014                         expanded = 1;
1015         }
1016         if (!mview->showHydrogens && an1 == 1)
1017                 return;
1018         if (!mview->showDummyAtoms && an1 == 0)
1019                 return;
1020         if (!mview->showExpandedAtoms && expanded)
1021                 return;
1022         if (ap != NULL && (ap->exflags & kAtomHiddenFlag))
1023                 return;
1024         dp = &(gElementParameters[an1]);
1025         if (dp == NULL)
1026                 return;
1027         if (selected) {
1028                 memcpy(rgba, sRedColor, sizeof(rgba));
1029         } else {
1030                 rgba[0] = dp->red / 65535.0;
1031                 rgba[1] = dp->green / 65535.0;
1032                 rgba[2] = dp->blue / 65535.0;
1033                 rgba[3] = 1.0;
1034         }
1035         if (expanded || periodicOffset != NULL) {
1036                 rgba[0] *= 0.5;
1037                 rgba[1] *= 0.5;
1038                 rgba[2] *= 0.5;
1039         }
1040         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
1041         if (periodicOffset != NULL)
1042                 VecInc(r1, *periodicOffset);
1043         p[0] = r1.x; p[1] = r1.y; p[2] = r1.z;
1044         if (mview->draggingMode == kMainViewDraggingSelectedAtoms && selected) {
1045                 p[0] += dragOffset->x;
1046                 p[1] += dragOffset->y;
1047                 p[2] += dragOffset->z;
1048         }
1049         rad = VDW_RADIUS(dp) * mview->atomRadius;
1050         if (mview->showEllipsoids) {
1051                 if (ap != NULL && ap->aniso != NULL) {
1052                         GLfloat elip[9];
1053                         Mat33 pmat2;
1054                         int i;
1055                         if (trp != NULL) {
1056                                 MatrixMul(pmat2, mview->mol->cell->rtr, ap->aniso->pmat);
1057                                 MatrixMul(pmat2, *((Mat33 *)trp), pmat2);
1058                                 MatrixMul(pmat2, mview->mol->cell->tr, pmat2);
1059                                 for (i = 0; i < 9; i++)
1060                                         elip[i] = pmat2[i] * mview->probabilityScale;
1061                         } else {
1062                                 for (i = 0; i < 9; i++)
1063                                         elip[i] = ap->aniso->pmat[i] * mview->probabilityScale;
1064                         }
1065                         drawEllipsoid(p, elip, elip+3, elip+6, mview->atomResolution * 3 / 2); /* Use higher resolution than spheres */
1066                 } else {
1067                         if (ap != NULL) {
1068                                 //  Recalculate radius from temperature factor
1069                                 rad = biso2radius(ap->tempFactor);
1070                                 rad *= mview->probabilityScale;
1071                         }
1072                         drawSphere(p, rad, mview->atomResolution);
1073                 }
1074         } else {
1075                 drawSphere(p, rad, mview->atomResolution);
1076         }
1077         pp[0] = p[0];
1078         pp[1] = p[1];
1079         pp[2] = p[2];
1080         if (MainView_convertObjectPositionToScreenPosition(mview, pp, p + 3)) {
1081         /*      fprintf(stderr, "atom %d: {%f, %f, %f}\n", i1, p[3], p[4], p[5]); */
1082                 float fp[3];
1083                 fp[0] = p[3]; fp[1] = p[4]; fp[2] = p[5];
1084                 MainViewCallback_drawLabel(mview, fp, label);
1085         }
1086 }
1087
1088 static void
1089 drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft, const Vector *dragOffset, const Vector *periodicOffset, int isAnchorBond)
1090 {
1091         const ElementPar *dp;
1092         int i, in;
1093         int an[2];
1094         char expanded[2];
1095         char anchor[2];
1096         Vector r[2];
1097         GLfloat p[6];
1098         GLfloat rgba[4];
1099         float rad_mul = 1.0;
1100         float alpha_mul = 1.0;
1101         int natoms = mview->mol->natoms;
1102         expanded[0] = expanded[1] = 0;
1103         anchor[0] = anchor[1] = 0;
1104
1105         for (i = 0; i < 2; i++) {
1106                 const Atom *ap;
1107                 in = (i == 0 ? i1 : i2);
1108                 if (in >= natoms && in < natoms + 2) {
1109                         if (mview->tempAtoms[in - natoms] >= 0) {
1110                                 ap = ATOM_AT_INDEX(mview->mol->atoms, mview->tempAtoms[in - natoms]);
1111                                 an[i] = ap->atomicNumber;
1112                                 r[i] = ap->r;
1113                         } else {
1114                                 ap = NULL;
1115                                 r[i] = mview->tempAtomPos[in - natoms];
1116                                 an[i] = 6;
1117                         }
1118                 } else {
1119                         ap = ATOM_AT_INDEX(mview->mol->atoms, in);
1120                         an[i] = ap->atomicNumber;
1121                         r[i] = ap->r;
1122                         if (SYMOP_ALIVE(ap->symop))
1123                                 expanded[i] = 1;
1124                         if (ap->anchor != NULL)
1125                                 anchor[i] = 1;
1126                 }
1127                 if (!mview->showHydrogens && an[i] == 1)
1128                         return;
1129                 if (!mview->showDummyAtoms && an[i] == 0)
1130                         return;
1131                 if (!mview->showExpandedAtoms && expanded[i])
1132                         return;
1133                 if (ap != NULL && (ap->exflags & kAtomHiddenFlag))
1134                         return;
1135         }
1136
1137         if (periodicOffset != NULL) {
1138                 VecInc(r[0], *periodicOffset);
1139                 VecInc(r[1], *periodicOffset);
1140         }
1141
1142         if (anchor[0] + anchor[1] == 2)
1143                 alpha_mul = 0.5;
1144         if (anchor[0] + anchor[1] == 1)
1145                 rad_mul = (isAnchorBond ? 0.3 : 0.6);
1146
1147         dp = &(gElementParameters[an[0]]);
1148         if (dp == NULL)
1149                 return;
1150         if (selected && selected2) {
1151                 memcpy(rgba, sRedColor, sizeof(rgba));
1152         } else {
1153                 rgba[0] = dp->red / 65535.0;
1154                 rgba[1] = dp->green / 65535.0;
1155                 rgba[2] = dp->blue / 65535.0;
1156                 rgba[3] = 1.0;
1157         }
1158         rgba[3] *= alpha_mul;
1159         if (expanded[0] || periodicOffset != NULL) {
1160                 rgba[0] *= 0.5;
1161                 rgba[1] *= 0.5;
1162                 rgba[2] *= 0.5;
1163         }
1164         if (mview->draggingMode == kMainViewDraggingSelectedAtoms) {
1165                 if (selected)
1166                         VecInc(r[0], *dragOffset);
1167                 if (selected2)
1168                         VecInc(r[1], *dragOffset);
1169         }
1170         p[0] = r[0].x; p[1] = r[0].y; p[2] = r[0].z;
1171         p[3] = (r[1].x + p[0]) * 0.5;
1172         p[4] = (r[1].y + p[1]) * 0.5;
1173         p[5] = (r[1].z + p[2]) * 0.5;
1174         if (draft) {
1175                 glColor3f(rgba[0], rgba[1], rgba[2]);
1176                 glBegin(GL_LINES);
1177                 glVertex3f(p[0], p[1], p[2]);
1178                 glVertex3f(p[3], p[4], p[5]);
1179                 glEnd();
1180         } else {
1181                 glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
1182                 drawCylinder(p, p + 3, mview->bondRadius * rad_mul, mview->bondResolution, 0);
1183         }
1184         dp = &(gElementParameters[an[1]]);
1185         if (dp == NULL)
1186                 return;
1187         if (!selected || !selected2) {
1188                 rgba[0] = dp->red / 65535.0;
1189                 rgba[1] = dp->green / 65535.0;
1190                 rgba[2] = dp->blue / 65535.0;
1191                 rgba[3] = 1.0;
1192         }
1193         rgba[3] *= alpha_mul;
1194         if (expanded[1] || periodicOffset != NULL) {
1195                 rgba[0] *= 0.5;
1196                 rgba[1] *= 0.5;
1197                 rgba[2] *= 0.5;
1198         }
1199         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
1200         p[0] = r[1].x; p[1] = r[1].y; p[2] = r[1].z;
1201         if (draft) {
1202                 glColor3f(rgba[0], rgba[1], rgba[2]);
1203                 glBegin(GL_LINES);
1204                 glVertex3f(p[0], p[1], p[2]);
1205                 glVertex3f(p[3], p[4], p[5]);
1206                 glEnd();
1207         } else {
1208                 glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
1209                 drawCylinder(p, p + 3, mview->bondRadius * rad_mul, 8, 0);
1210         }
1211 }
1212
1213 /*  Calculate drag offset during moving the selection.  */
1214 static void
1215 calcDragOffset(MainView *mview, Vector *outVector)
1216 {
1217         double p[6];
1218         if (mview->draggingMode == kMainViewDraggingSelectedAtoms
1219         && MainView_convertScreenPositionToObjectPosition(mview, mview->dragStartPos, p)
1220         && MainView_convertScreenPositionToObjectPosition(mview, mview->dragEndPos, p + 3)) {
1221                 outVector->x = p[3] - p[0];
1222                 outVector->y = p[4] - p[1];
1223                 outVector->z = p[5] - p[2];
1224         } else {
1225                 outVector->x = outVector->y = outVector->z = 0;
1226         }
1227
1228 }
1229
1230 static void
1231 drawSurface(MainView *mview)
1232 {
1233         int i, sn, k;
1234         GLfloat rgba[4];
1235         MCube *mc;
1236         if (mview->mol == NULL || mview->mol->mcube == NULL)
1237                 return;
1238         mc = mview->mol->mcube;
1239         for (sn = 0; sn <= 1; sn++) {
1240                 if (mc->c[sn].ntriangles == 0)
1241                         continue;
1242                 for (i = 0; i < 4; i++)
1243                         rgba[i] = mc->c[sn].rgba[i];
1244                 k = (sn == 0 ? -1 : 1);
1245                 glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
1246                 glBegin(GL_TRIANGLES);
1247                 for (i = 0; mc->c[sn].triangles[i] >= 0; i++) {
1248                         MCubePoint *mcp = mc->c[sn].cubepoints + mc->c[sn].triangles[i];
1249                         glNormal3f(mcp->grad[0] * k, mcp->grad[1] * k, mcp->grad[2] * k);
1250                         glVertex3f(mcp->pos[0], mcp->pos[1], mcp->pos[2]);
1251                 }
1252                 glEnd();
1253         }
1254 }
1255
1256 static void
1257 drawModel(MainView *mview)
1258 {
1259         Molecule *mol;
1260     int i, natoms, nbonds;
1261         int amin, amax, bmin, bmax, cmin, cmax, da, db, dc;
1262         Byte original;
1263         Double atomRadius, bondRadius;
1264         Vector periodicOffset;
1265         Vector *axes;
1266         int selected, selected2;
1267         char *selectFlags;
1268         Atom *ap;
1269         int draft = mview->lineMode;
1270 /*    static Double gray[] = {0.8, 0.8, 0.8, 1}; */
1271
1272         atomRadius = mview->atomRadius;
1273         bondRadius = mview->bondRadius;
1274         mol = mview->mol;
1275         natoms = mol->natoms;
1276 /*    if (natoms == 0) {
1277         return;
1278     }
1279 */
1280         if (mview->draggingMode == kMainViewSelectingRegion)
1281                 selectFlags = temporarySelection(mview, mview->modifierFlags, 0, 0);
1282         else selectFlags = NULL;
1283
1284         if (mview->draggingMode == kMainViewDraggingSelectedAtoms)
1285                 calcDragOffset(mview, &mview->dragOffset);
1286         else mview->dragOffset.x = mview->dragOffset.y = mview->dragOffset.z = 0;
1287
1288         if (mview->showGraphite > 0 && mview->showGraphiteFlag) {
1289                 drawGraphite(mview);
1290         }
1291
1292         amin = amax = bmin = bmax = cmin = cmax = 0;
1293         if (mview->showExpandedAtoms && mol->cell != NULL) {
1294                 if (mol->cell->flags[0] && mview->showPeriodicImageFlag) {
1295                         amin = mview->showPeriodicImage[0];
1296                         amax = mview->showPeriodicImage[1];
1297                 }
1298                 if (mol->cell->flags[1] && mview->showPeriodicImageFlag) {
1299                         bmin = mview->showPeriodicImage[2];
1300                         bmax = mview->showPeriodicImage[3];
1301                 }
1302                 if (mol->cell->flags[2] && mview->showPeriodicImageFlag) {
1303                         cmin = mview->showPeriodicImage[4];
1304                         cmax = mview->showPeriodicImage[5];
1305                 }
1306                 axes = mol->cell->axes;
1307         } else {
1308                 axes = NULL;
1309         }
1310
1311         if (draft == 0) {
1312                 for (da = amin; da <= amax; da++) {
1313                         for (db = bmin; db <= bmax; db++) {
1314                                 for (dc = cmin; dc <= cmax; dc++) {
1315                                         original = (da == 0 && db == 0 && dc == 0);
1316                                         if (!original) {
1317                                                 VecScale(periodicOffset, axes[0], da);
1318                                                 VecScaleInc(periodicOffset, axes[1], db);
1319                                                 VecScaleInc(periodicOffset, axes[2], dc);
1320                                         }
1321                                         for (i = 0, ap = mview->mol->atoms; i < natoms; i++, ap = ATOM_NEXT(ap)) {
1322                                                 if (mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
1323                                                         /*  Mouse event is detected  */
1324                                                         draft = 1;
1325                                                         goto skip;
1326                                                 }
1327                                                 if (mview->draggingMode == kMainViewCreatingBond && (i == mview->tempAtoms[0] || i == mview->tempAtoms[1] || i >= natoms))
1328                                                         selected = 1;  /*  extra atoms  */
1329                                                 else if (selectFlags != NULL)
1330                                                         selected = selectFlags[i];
1331                                                 else
1332                                                         selected = MoleculeIsAtomSelected(mview->mol, i);
1333                                                 drawAtom(mview, i, selected, &mview->dragOffset, (original ? NULL : &periodicOffset));
1334                                                 if (ap->anchor != NULL) {
1335                                                         /*  Draw anchor bonds  */
1336                                                         Int j, k;
1337                                                         Int *cp = AtomConnectData(&ap->anchor->connect);
1338                                                         for (j = 0; j < ap->anchor->connect.count; j++) {
1339                                                                 k = cp[j];
1340                                                                 if (selectFlags != NULL)
1341                                                                         selected2 = selectFlags[k];
1342                                                                 else
1343                                                                         selected2 = MoleculeIsAtomSelected(mview->mol, k);
1344                                                                 drawBond(mview, i, k, selected, selected2, draft, &mview->dragOffset, (original ? NULL : &periodicOffset), 1);
1345                                                         }
1346                                                 }
1347                                         }
1348                                 }
1349
1350                                 if (draft == 0) {
1351                                         if (original) {
1352                                                 /*  Extra atoms  */
1353                                                 drawAtom(mview, natoms, 1, &mview->dragOffset, NULL);
1354                                                 drawAtom(mview, natoms + 1, 1, &mview->dragOffset, NULL);
1355                                         }
1356                                         /*  Expanded atoms  */
1357                                 /*      if (mview->showExpandedAtoms) {
1358                                                 for (i = 0; i < mview->mol->nexatoms; i++) {
1359                                                         if (mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
1360                                                                 //  Mouse event is detected
1361                                                                 draft = 1;
1362                                                                 break;
1363                                                         //      goto cleanup;
1364                                                         }
1365                                                         drawAtom(mview, -i-1, (selectFlags == NULL ? 0 : selectFlags[i]), &dragOffset, (original ? NULL : &periodicOffset));
1366                                                 }
1367                                         } */
1368                                 }
1369                         }
1370                 }
1371         }
1372
1373 skip:
1374         nbonds = mol->nbonds;
1375         if (draft)
1376                 glDisable(GL_LIGHTING);
1377         for (da = amin; da <= amax; da++) {
1378                 for (db = bmin; db <= bmax; db++) {
1379                         for (dc = cmin; dc <= cmax; dc++) {
1380                                 original = (da == 0 && db == 0 && dc == 0);
1381                                 if (!original) {
1382                                         VecScale(periodicOffset, axes[0], da);
1383                                         VecScaleInc(periodicOffset, axes[1], db);
1384                                         VecScaleInc(periodicOffset, axes[2], dc);
1385                                 }
1386
1387                                 for (i = 0; i < nbonds; i++) {
1388                                         int n1, n2;
1389                                         if (draft == 0 && mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
1390                                                 /*  Mouse event is detected  */
1391                                                 draft = 1;
1392                                                 glDisable(GL_LIGHTING);
1393                                         /*      goto cleanup;  */
1394                                         }
1395                                         n1 = mview->mol->bonds[i * 2];
1396                                         n2 = mview->mol->bonds[i * 2 + 1];
1397                                         if (selectFlags == NULL) {
1398                                                 selected = MoleculeIsAtomSelected(mview->mol, n1);
1399                                                 selected2 = MoleculeIsAtomSelected(mview->mol, n2);
1400                                         } else {
1401                                                 selected = selectFlags[n1];
1402                                                 selected2 = selectFlags[n2];
1403                                         }
1404                                         drawBond(mview, n1, n2, selected, selected2, draft, &mview->dragOffset, (original ? NULL : &periodicOffset), 0);
1405                                 }
1406
1407                                 /*  Extra bond  */
1408                                 if (original && mview->draggingMode == kMainViewCreatingBond) {
1409                                         drawBond(mview, natoms, natoms + 1, 1, 1, draft, &mview->dragOffset, NULL, 0);
1410                                 }
1411
1412                                 /*  Expanded bonds  */
1413                         /*      for (i = 0; i < mview->mol->nexbonds; i++) {
1414                                         int n1, n2;
1415                                         if (draft == 0 && mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
1416                                                 //  Mouse event is detected
1417                                                 draft = 1;
1418                                                 glDisable(GL_LIGHTING);
1419                                         //      goto cleanup;
1420                                         }
1421                                         n1 = mview->mol->exbonds[i * 2];
1422                                         n2 = mview->mol->exbonds[i * 2 + 1];
1423                                         if (n1 < 0)
1424                                                 n1 = -n1 - 1;
1425                                         if (n2 < 0)
1426                                                 n2 = -n2 - 1;
1427                                         if (selectFlags == NULL) {
1428                                                 selected = MoleculeIsAtomSelected(mview->mol, n1);
1429                                                 selected2 = MoleculeIsAtomSelected(mview->mol, n2);
1430                                         } else {
1431                                                 selected = selectFlags[n1];
1432                                                 selected2 = selectFlags[n2];
1433                                         }
1434                                         drawBond(mview, mview->mol->exbonds[i * 2], mview->mol->exbonds[i * 2 + 1], selected, selected2, draft, &dragOffset, (original ? NULL : &periodicOffset));
1435                                 } */
1436                         }
1437                 }
1438         }
1439
1440 /*  cleanup: */
1441         if (draft)
1442                 enableLighting();
1443         if (selectFlags != NULL)
1444                 free(selectFlags);
1445 }
1446
1447 static void
1448 drawGraphics(MainView *mview)
1449 {
1450         int i, j;
1451         MainViewGraphic *g;
1452         for (i = 0; i < mview->ngraphics; i++) {
1453                 g = &mview->graphics[i];
1454                 if (g->visible == 0)
1455                         continue;
1456                 glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, g->rgba);
1457                 switch (g->kind) {
1458                         case kMainViewGraphicLine:
1459                                 glDisable(GL_LIGHTING);
1460                                 glColor4fv(g->rgba);
1461                                 glBegin(GL_LINE_STRIP);
1462                                 for (j = 0; j < g->npoints; j++) {
1463                                         if (g->points[j * 3] >= kInvalidFloat)
1464                                                 break;
1465                                         glVertex3fv(&g->points[j * 3]);
1466                                 }
1467                                 glEnd();
1468                                 enableLighting();
1469                                 break;
1470                         case kMainViewGraphicPoly: {
1471                                 Vector v0, v1, v2, v3;
1472                                 glBegin(GL_TRIANGLE_FAN);
1473                                 v1.x = g->points[0] - g->points[g->npoints - 3];
1474                                 v1.y = g->points[1] - g->points[g->npoints - 2];
1475                                 v1.z = g->points[2] - g->points[g->npoints - 1];
1476                                 v0 = v1;
1477                                 for (j = 0; j < g->npoints; j++) {
1478                                         v2.x = g->points[j * 3 + 3] - g->points[j * 3];
1479                                         v2.y = g->points[j * 3 + 4] - g->points[j * 3 + 1];
1480                                         v2.z = g->points[j * 3 + 5] - g->points[j * 3 + 2];
1481                                         VecCross(v3, v1, v2);
1482                                         if (NormalizeVec(&v3, &v3) == 0)
1483                                                 glNormal3f(v3.x, v3.y, v3.z);
1484                                         glVertex3fv(&g->points[j * 3]);
1485                                         v1 = v2;
1486                                 }
1487                                 if (g->closed) {
1488                                         VecCross(v3, v1, v0);
1489                                         if (NormalizeVec(&v3, &v3) == 0)
1490                                                 glNormal3f(v3.x, v3.y, v3.z);
1491                                         glVertex3fv(g->points);
1492                                 }
1493                                 glEnd();
1494                                 break;
1495                         }
1496                         case kMainViewGraphicCylinder:
1497                                 drawCylinder(g->points, g->points + 3, g->points[6], 15, g->closed);
1498                                 break;
1499                         case kMainViewGraphicCone:
1500                                 drawCone(g->points, g->points + 3, g->points[6], 15, g->closed);
1501                                 break;
1502                         case kMainViewGraphicEllipsoid:
1503                                 drawEllipsoid(g->points, g->points + 3, g->points + 6, g->points + 9, 8);
1504                                 break;
1505                 }
1506         }
1507 }
1508
1509 static void
1510 drawUnitCell(MainView *mview)
1511 {
1512         GLfloat a[3], b[3], c[3], ab[3], bc[3], ca[3], abc[3];
1513         XtalCell *cp;
1514         GLfloat origin[3];
1515         if (!mview->showUnitCell || (cp = mview->mol->cell) == NULL)
1516                 return;
1517         origin[0] = cp->origin.x;
1518         origin[1] = cp->origin.y;
1519         origin[2] = cp->origin.z;
1520         a[0] = cp->axes[0].x + origin[0];
1521         a[1] = cp->axes[0].y + origin[1];
1522         a[2] = cp->axes[0].z + origin[2];
1523         b[0] = cp->axes[1].x + origin[0];
1524         b[1] = cp->axes[1].y + origin[1];
1525         b[2] = cp->axes[1].z + origin[2];
1526         c[0] = cp->axes[2].x + origin[0];
1527         c[1] = cp->axes[2].y + origin[1];
1528         c[2] = cp->axes[2].z + origin[2];
1529
1530         ab[0] = a[0] + b[0]; ab[1] = a[1] + b[1]; ab[2] = a[2] + b[2];
1531         bc[0] = b[0] + c[0]; bc[1] = b[1] + c[1]; bc[2] = b[2] + c[2];
1532         ca[0] = c[0] + a[0]; ca[1] = c[1] + a[1]; ca[2] = c[2] + a[2];
1533         abc[0] = a[0] + bc[0]; abc[1] = a[1] + bc[1]; abc[2] = a[2] + bc[2];
1534
1535         glDisable(GL_LIGHTING);
1536         glColor3f(0.75, 0.2, 0.0);
1537         glBegin(GL_LINES);
1538         glVertex3fv(origin);
1539         glVertex3fv(a);
1540         glVertex3fv(origin);
1541         glVertex3fv(b);
1542         glVertex3fv(origin);
1543         glVertex3fv(c);
1544         glVertex3fv(a);
1545         glVertex3fv(ab);
1546         glVertex3fv(a);
1547         glVertex3fv(ca);
1548         glVertex3fv(b);
1549         glVertex3fv(ab);
1550         glVertex3fv(b);
1551         glVertex3fv(bc);
1552         glVertex3fv(c);
1553         glVertex3fv(ca);
1554         glVertex3fv(c);
1555         glVertex3fv(bc);
1556         glVertex3fv(ab);
1557         glVertex3fv(abc);
1558         glVertex3fv(bc);
1559         glVertex3fv(abc);
1560         glVertex3fv(ca);
1561         glVertex3fv(abc);
1562         glEnd();
1563
1564         enableLighting();
1565 }
1566
1567 static void
1568 drawRotationCenter(MainView *mview)
1569 {
1570         GLfloat ps[3], pe[3], col[4];
1571         double fd[2];  /*  Fovy and distance  */
1572         double tr[3];  /*  Translation  */
1573         double r, rr;
1574         if (mview == NULL || !mview->showRotationCenter)
1575                 return;
1576         TrackballGetTranslate(mview->track, tr);
1577         TrackballGetPerspective(mview->track, fd);
1578         tr[0] *= -mview->dimension;
1579         tr[1] *= -mview->dimension;
1580         tr[2] *= -mview->dimension;
1581         r = fd[1] * mview->dimension * tan(fd[0] * 0.5 * kDeg2Rad) * 0.1;
1582         rr = r * 0.1;
1583         ps[0] = tr[0];
1584         ps[1] = tr[1];
1585         ps[2] = tr[2];
1586         col[0] = col[1] = col[2] = 0.5; col[3] = 1.0;
1587         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
1588         drawSphere(ps, rr, 8);
1589         ps[0] = tr[0] - r;
1590         pe[0] = tr[0] + r;
1591         pe[1] = tr[1];
1592         pe[2] = tr[2];
1593         col[0] = 1.0; col[1] = col[2] = 0.0;
1594         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
1595         drawSphere(ps, rr, 8);
1596         drawSphere(pe, rr, 8);
1597         drawCylinder(ps, pe, rr, 8, 0);
1598         ps[0] = tr[0];
1599         ps[1] = tr[1] - r;
1600         pe[0] = tr[0];
1601         pe[1] = tr[1] + r;
1602         col[1] = 1.0; col[0] = col[2] = 0.0;
1603         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
1604         drawSphere(ps, rr, 8);
1605         drawSphere(pe, rr, 8);
1606         drawCylinder(ps, pe, rr, 8, 0);
1607         ps[1] = tr[1];
1608         ps[2] = tr[2] - r;
1609         pe[1] = tr[1];
1610         pe[2] = tr[2] + r;
1611         col[2] = 1.0; col[0] = col[1] = 0.0;
1612         glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
1613         drawSphere(ps, rr, 8);
1614         drawSphere(pe, rr, 8);
1615         drawCylinder(ps, pe, rr, 8, 0);
1616 }
1617
1618 static int
1619 compareLabelByDepth(const void *ap, const void *bp)
1620 {
1621         Double dz = ((const LabelRecord *)bp)->pos.z - ((const LabelRecord *)ap)->pos.z;
1622         if (dz > 0)
1623                 return 1;
1624         else if (dz < 0)
1625                 return -1;
1626         else return 0;
1627 }
1628
1629 static void
1630 drawLabels(MainView *mview)
1631 {
1632 /*      Transform *trp; */
1633         Atom *ap;
1634         LabelRecord *lp;
1635         int i, nlabels;
1636
1637         if (mview->nlabels == 0)
1638                 return;
1639
1640 /*      mview->sortedLabels = (LabelRecord **)calloc(sizeof(LabelRecord *), mview->nlabels);
1641         if (mview->sortedLabels == NULL)
1642                 return; */
1643
1644 /*      if (mview->mol->is_xtal_coord)
1645                 trp = &(mview->mol->cell->tr);
1646         else trp = NULL; */
1647
1648         /*  Get the screen coordinates of the labels  */
1649         nlabels = 0;
1650         for (i = 0; i < mview->mol->natoms; i++) {
1651                 float scrp[3], f[3];
1652                 ap = ATOM_AT_INDEX(mview->mol->atoms, i);
1653                 if (ap->exflags & kAtomHiddenFlag)
1654                         continue;
1655                 if (!mview->showHydrogens && ap->atomicNumber == 1)
1656                         continue;
1657                 if (!mview->showDummyAtoms && ap->atomicNumber == 0)
1658                         continue;
1659                 if (ap->labelid <= 0 || ap->labelid > mview->nlabels)
1660                         continue;
1661                 lp = mview->labels + (ap->labelid - 1);
1662                 MainView_screenCenterPointOfAtom(mview, lp->idx1, scrp);
1663         /*      if (lp->idx2 >= 0) {
1664                         Atom *bp;
1665                         r1 = ap->r;
1666                         if (trp != NULL)
1667                                 TransformVec(&r1, *trp, &r1);
1668                         bp = ATOM_AT_INDEX(mview->mol->atoms, lp->idx2);
1669                         if (bp->exflags & kAtomHiddenFlag)
1670                                 continue;
1671                         r2 = bp->r;
1672                         if (trp != NULL)
1673                                 TransformVec(&r2, *trp, &r2);
1674                         r1.x = (r1.x + r2.x) * 0.5;
1675                         r1.y = (r1.y + r2.y) * 0.5;
1676                         r1.z = (r1.z + r2.z) * 0.5;
1677                         objp[0] = r1.x;
1678                         objp[1] = r1.y;
1679                         objp[2] = r1.z;
1680                         MainView_convertObjectPositionToScreenPosition(mview, objp, scrp);
1681                 } */
1682                 lp->pos.x = scrp[0];
1683                 lp->pos.y = scrp[1];
1684                 lp->pos.z = scrp[2];
1685                 MainViewCallback_labelSize(lp->label, f);
1686                 f[0] = floor(lp->pos.x - f[0] * 0.5);
1687                 f[1] = -floor(lp->pos.y + f[1] * 0.5);
1688                 f[2] = lp->pos.z;
1689         /*      fprintf(stderr, "label position (%d) = {%f, %f, %f}\n", i, f[0], f[1], f[2]); */
1690                 MainViewCallback_drawLabelAtPoint(lp->label, f);
1691         /*      mview->sortedLabels[nlabels++] = lp;
1692                 if (nlabels >= mview->nlabels)
1693                         break; */
1694         }
1695
1696 /*      //  Sort labels by z coordinates (descending order)
1697         qsort(mview->sortedLabels, nlabels, sizeof(LabelRecord *), compareLabelByDepth);
1698
1699         //  Draw labels
1700         for (i = 0; i < nlabels; i++) {
1701                 LabelRecord *lp = mview->sortedLabels[i];
1702                 Double f[3];
1703                 MainViewCallback_labelSize(lp->label, f);
1704                 f[0] = floor(lp->pos.x - f[0] * 0.5);
1705                 f[1] = -floor(lp->pos.y + f[1] * 0.5);
1706                 f[2] = lp->pos.z;
1707                 MainViewCallback_drawLabelAtPoint(lp->label, f);
1708         }
1709 */
1710 }
1711
1712 void
1713 MainView_drawModel(MainView *mview)
1714 {
1715     double w[4], dimension, distance;
1716         float frame[4], width, height, scale;
1717         GLdouble *pp;
1718         Transform mtr;
1719
1720 /*    int i;  */
1721
1722         if (mview == NULL)
1723                 return;
1724
1725     if (!mview->isInitialized) {
1726         MainView_initializeOpenGL();
1727                 mview->isInitialized = 1;
1728     }
1729
1730         dimension = mview->dimension;
1731
1732         if (mview->offline_scale == 0.0)
1733                 scale = 1.0;
1734         else
1735                 scale = mview->offline_scale;
1736
1737         MainViewCallback_frame(mview, frame);
1738         width = (frame[2] - frame[0]) * scale;
1739         height = (frame[3] - frame[1]) * scale;
1740
1741     glViewport(0, 0, width, height);
1742
1743         /*  Clear the buffer  */
1744     glClearColor(mview->background_color[0], mview->background_color[1], mview->background_color[2], mview->background_color[3]);
1745     glClear(GL_COLOR_BUFFER_BIT |
1746             GL_DEPTH_BUFFER_BIT);
1747
1748         if (mview->mol == NULL)
1749                 return;
1750
1751         /*  Set up the projection  */
1752     glMatrixMode(GL_PROJECTION);
1753     glLoadIdentity();
1754         TrackballGetPerspective(mview->track, w);
1755     distance = w[1] * dimension;
1756     mview->perspective_vector[0] = w[0];
1757     mview->perspective_vector[1] = width / height;
1758     mview->perspective_vector[2] = dimension;
1759     mview->perspective_vector[3] = distance + 200.0 * dimension;
1760     gluPerspective(mview->perspective_vector[0], mview->perspective_vector[1], mview->perspective_vector[2], mview->perspective_vector[3]);
1761
1762     /*  Set up the model view  */
1763     glMatrixMode(GL_MODELVIEW);
1764     glLoadIdentity();
1765     glTranslatef(0.0, 0.0, -distance);
1766         TrackballGetRotate(mview->track, w);
1767     glRotatef(w[0], w[1], w[2], w[3]);
1768         TrackballGetTranslate(mview->track, w);
1769         w[0] *= dimension;
1770         w[1] *= dimension;
1771         w[2] *= dimension;
1772     glTranslatef(w[0], w[1], w[2]);
1773         mview->lookat.x = -w[0];
1774         mview->lookat.y = -w[1];
1775         mview->lookat.z = -w[2];
1776
1777         MainViewCallback_clearLabels(mview);
1778     drawModel(mview);
1779         drawSurface(mview);
1780         drawUnitCell(mview);
1781         drawRotationCenter(mview);
1782         drawGraphics(mview);
1783
1784         /*  Get important matrices and vectors  */
1785     glGetDoublev(GL_MODELVIEW_MATRIX, mview->modelview_matrix);
1786     glGetDoublev(GL_PROJECTION_MATRIX, mview->projection_matrix);
1787         pp = mview->modelview_matrix;
1788         mtr[0] = pp[0]; mtr[1] = pp[1]; mtr[2] = pp[2];
1789         mtr[3] = pp[4]; mtr[4] = pp[5]; mtr[5] = pp[6];
1790         mtr[6] = pp[8]; mtr[7] = pp[9]; mtr[8] = pp[10];
1791         mtr[9] = pp[12]; mtr[10] = pp[13]; mtr[11] = pp[14];
1792         TransformInvert(mtr, mtr);
1793         mview->camera.x = mtr[9]; mview->camera.y = mtr[10]; mview->camera.z = mtr[11];
1794         mview->lookto.x = mtr[6]; mview->lookto.y = mtr[7]; mview->lookto.z = mtr[8];
1795         mview->up.x = mtr[3]; mview->up.y = mtr[4]; mview->up.z = mtr[5];
1796
1797         /*  Draw labels  */
1798         glDisable(GL_LIGHTING);
1799 //      glDisable (GL_DEPTH_TEST);
1800 //      glEnable (GL_BLEND);
1801 //      glBlendFunc (GL_ONE, GL_ONE_MINUS_SRC_ALPHA);
1802 #if __WXMAC__
1803         glEnable (GL_TEXTURE_RECTANGLE_EXT);
1804 #endif
1805         glMatrixMode (GL_PROJECTION);
1806         glLoadIdentity ();
1807         glMatrixMode (GL_MODELVIEW);
1808         glLoadIdentity ();
1809         glOrtho(0, width, 0, -height, 0.0, -1.0);  /*  non-flipped view  */
1810         drawLabels(mview);
1811 #if __WXMAC__
1812         glDisable (GL_TEXTURE_RECTANGLE_EXT);
1813 #endif
1814 //      glDisable(GL_BLEND);
1815 //      glEnable(GL_DEPTH_TEST);
1816         enableLighting();
1817
1818     if (mview->draggingMode == kMainViewSelectingRegion) {
1819                 /*  Draw selection rectangle  */
1820                 glDisable(GL_LIGHTING);
1821                 glDisable(GL_DEPTH_TEST);
1822         glMatrixMode(GL_MODELVIEW);
1823         glLoadIdentity();
1824         glMatrixMode(GL_PROJECTION);
1825         glLoadIdentity();
1826         glOrtho(0, width, 0, height, -1.0, 1.0);
1827         glColor3f(1.0, 1.0, 0.0);
1828         glBegin(GL_LINE_STRIP);
1829                 glVertex2f(mview->dragStartPos[0], mview->dragStartPos[1]);
1830                 glVertex2f(mview->dragStartPos[0], mview->dragEndPos[1]);
1831                 glVertex2f(mview->dragEndPos[0], mview->dragEndPos[1]);
1832                 glVertex2f(mview->dragEndPos[0], mview->dragStartPos[1]);
1833                 glVertex2f(mview->dragStartPos[0], mview->dragStartPos[1]);
1834         glEnd();
1835                 glEnable(GL_DEPTH_TEST);
1836                 enableLighting();
1837     }
1838
1839     glFinish();
1840
1841 }
1842
1843 #pragma mark ====== Labels ======
1844
1845 void
1846 MainView_attachLabelToAtom(MainView *mview, int index)
1847 {
1848         LabelRecord rec;
1849         char buf[24], *p;
1850         Atom *ap;
1851 /*      const ElementPar *dp; */
1852         static float foreColor[] = {1, 1, 0, 1};
1853         static float backColor[] = {1, 1, 1, 0};
1854         ap = ATOM_AT_INDEX(mview->mol->atoms, index);
1855         if (ap->resSeq == 0)
1856                 snprintf(buf, sizeof buf, "%-.4s", ap->aname);
1857         else
1858                 snprintf(buf, sizeof buf, "%d:%-.4s", ap->resSeq, ap->aname);
1859         for (p = buf; *p; p++) {
1860                 if (isspace(*p)) {
1861                         *p = 0;
1862                         break;
1863                 }
1864         }
1865 /*      dp = &(gBuiltinParameters->atomPars[ap->atomicNumber]);
1866         foreColor[0] = 1.0 - dp->r;
1867         foreColor[1] = 1.0 - dp->g;
1868         foreColor[2] = 1.0 - dp->b; */
1869         rec.label = MainViewCallback_newLabel(mview, buf, 14, foreColor, backColor);
1870         rec.idx1 = index;
1871         rec.idx2 = kInvalidIndex;
1872         rec.labelid = mview->nlabels + 1;
1873         ap->labelid = rec.labelid;
1874         AssignArray(&mview->labels, &mview->nlabels, sizeof(LabelRecord), mview->nlabels, &rec);
1875 }
1876
1877 void
1878 MainView_detachLabelFromAtom(MainView *mview, int index)
1879 {
1880         Atom *ap;
1881         if (index >= 0 && index < mview->mol->natoms) {
1882                 ap = ATOM_AT_INDEX(mview->mol->atoms, index);
1883                 if (ap->labelid > 0 && ap->labelid <= mview->nlabels) {
1884                         mview->labels[ap->labelid - 1].idx1 = kInvalidIndex;
1885                 }
1886                 ap->labelid = 0;
1887         }
1888 }
1889
1890 void
1891 MainView_purgeUnusedLabels(MainView *mview)
1892 {
1893         int i, n;
1894         Atom *ap;
1895         Int *tempid;
1896         if (mview == NULL || mview->nlabels == 0 || mview->mol->natoms == 0)
1897                 return;
1898         tempid = (Int *)calloc(sizeof(Int), mview->nlabels);
1899         if (tempid == NULL)
1900                 return;
1901
1902         /*  Mark the labels in use  */
1903         for (i = 0, ap = mview->mol->atoms; i < mview->mol->natoms; i++, ap = ATOM_NEXT(ap)) {
1904                 if (ap->labelid > 0 && ap->labelid <= mview->nlabels)
1905                         tempid[ap->labelid - 1] = 1;
1906         }
1907
1908         /*  Release the labels not in use, and assign new IDs  */
1909         n = 1;
1910         for (i = 0; i < mview->nlabels; i++) {
1911                 if (tempid[i] == 0) {
1912                         MainViewCallback_releaseLabel(mview->labels[i].label);
1913                         mview->labels[i].label = NULL;
1914                 } else {
1915                         mview->labels[i].labelid = tempid[i] = n++;
1916                 }
1917         }
1918
1919         /*  Purge the unused entries  */
1920         for (i = mview->nlabels - 1; i >= 0; i--) {
1921                 if (tempid[i] == 0) {
1922                         memmove(mview->labels + i + 1, mview->labels + i, (mview->nlabels - 1 - i) * sizeof(LabelRecord));
1923                         mview->nlabels--;
1924                 }
1925         }
1926         if (mview->nlabels == 0) {
1927                 free(mview->labels);
1928                 mview->labels = NULL;
1929         }
1930
1931         /*  Renumber  */
1932         for (i = 0, ap = mview->mol->atoms; i < mview->mol->natoms; i++, ap = ATOM_NEXT(ap)) {
1933                 if (ap->labelid > 0 && ap->labelid <= mview->nlabels) {
1934                         ap->labelid = tempid[ap->labelid - 1];
1935                 } else {
1936                         ap->labelid = 0;
1937                 }
1938         }
1939 }
1940
1941 #pragma mark ====== Mode ======
1942
1943 void
1944 MainView_setMode(MainView *mview, int mode)
1945 {
1946         if (mview != NULL)
1947                 mview->mode = mode;
1948 }
1949
1950 int
1951 MainView_getMode(const MainView *mview)
1952 {
1953         if (mview != NULL)
1954                 return mview->mode;
1955         else return 0;
1956 }
1957
1958 #pragma mark ====== Mouse operations ======
1959
1960 static void
1961 mousePosToTrackballPos(MainView *mview, const float *mousePos, float *pos)
1962 {
1963         float frame[4], width, height, radius;
1964         NULL_CHECK(mview, "mousePosToTrackballPos");
1965         MainViewCallback_frame(mview, frame);
1966         width = frame[2] - frame[0];
1967         height = frame[3] - frame[1];
1968         radius = (width > height ? height * 0.5 : width * 0.5);
1969         pos[0] = (mousePos[0] - frame[0] - width * 0.5) / radius;
1970         pos[1] = (mousePos[1] - frame[1] - height * 0.5) / radius;
1971 }
1972
1973 static void
1974 showAtomsInInfoText(MainView *mview, int n1, int n2)
1975 {
1976         char buf[64];
1977         if (n1 >= 0) {
1978                 MoleculeGetAtomName(mview->mol, n1, buf, sizeof buf - 1);
1979                 if (n2 >= 0) {
1980                         int nn = strlen(buf);
1981                         buf[nn++] = '-';
1982                         MoleculeGetAtomName(mview->mol, n2, buf + nn, sizeof buf - nn);
1983                 }
1984                 MainViewCallback_drawInfoText(mview, buf);
1985         } else MainViewCallback_drawInfoText(mview, NULL);
1986 }
1987
1988 void
1989 MainView_mouseDown(MainView *mview, const float *mousePos, int flags)
1990 {
1991         double p[3];
1992         float fp[3];
1993         float screenPos[3];
1994         double objectPos[3];
1995         int n1, n2, found;
1996         Atom *ap;
1997         mview->dragStartPos[0] = mousePos[0];
1998         mview->dragStartPos[1] = mousePos[1];
1999         mview->dragStartPos[2] = 0.5;
2000         found = MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 0, 0);
2001         mview->clickedAtoms[0] = n1;
2002         mview->clickedAtoms[1] = n2;
2003         if (found) {
2004                 /*  Estimate the screen z-coordinate of the mouse position  */
2005                 Vector r1, r2;
2006                 ap = ATOM_AT_INDEX(mview->mol->atoms, n1);
2007                 r1 = ap->r;
2008         /*      if (mview->mol->is_xtal_coord)
2009                                 TransformVec(&r1, mview->mol->cell->tr, &r1); */
2010                 objectPos[0] = r1.x;
2011                 objectPos[1] = r1.y;
2012                 objectPos[2] = r1.z;
2013                 if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, screenPos)) {
2014                         if (n2 >= 0) {
2015                                 ap = ATOM_AT_INDEX(mview->mol->atoms, n2);
2016                                 r2 = ap->r;
2017                         /*      if (mview->mol->is_xtal_coord)
2018                                                 TransformVec(&r2, mview->mol->cell->tr, &r2); */
2019                                 objectPos[0] = r2.x;
2020                                 objectPos[1] = r2.y;
2021                                 objectPos[2] = r2.z;
2022                                 if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, fp)) {
2023                                         Double w;
2024                                         r1.x = fp[0] - screenPos[0];
2025                                         r1.y = fp[1] - screenPos[1];
2026                                         r1.z = fp[2] - screenPos[2];
2027                                         NormalizeVec(&r1, &r1);
2028                                         r2.x = mousePos[0] - screenPos[0];
2029                                         r2.y = mousePos[1] - screenPos[1];
2030                                         r2.z = 0.5;
2031                                         w = VecDot(r1, r2);
2032                                         VecScale(r2, r1, w);
2033                                         screenPos[2] += r2.z;
2034                                 }
2035                         }
2036                 }
2037                 mview->dragStartPos[2] = screenPos[2];
2038         } else {
2039                 /*  Set the 'depth' value of the molecule center (i.e. trackball->translate * dimension)
2040                     as the z-coordinate of the drag start position  */
2041                 TrackballGetTranslate(mview->track, p);
2042                 objectPos[0] = -mview->dimension * p[0];
2043                 objectPos[1] = -mview->dimension * p[1];
2044                 objectPos[2] = -mview->dimension * p[2];
2045                 if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, screenPos))
2046                         mview->dragStartPos[2] = screenPos[2];
2047         }
2048
2049         mview->isDragging = 0;
2050
2051         switch (mview->mode) {
2052                 case kTrackballRotateMode:
2053                 case kTrackballTranslateMode:
2054                 case kTrackballScaleMode:
2055                         if (!found) {
2056                                 mousePosToTrackballPos(mview, mousePos, fp);
2057                                 TrackballStartDragging(mview->track, fp, mview->mode);
2058                                 mview->draggingMode = kMainViewMovingTrackball;
2059                                 break;
2060                         } else {
2061                                 /*  No 'break' intentional; drop to next case  */
2062                         }
2063                 case kTrackballSelectionMode: {
2064                         /*  Select or move around */
2065                         if (found) {
2066                                 if (flags & kShiftKeyMask) {
2067                                         /*  Shift-key pressed: toggle selection  */
2068                                         MoleculeToggleSelectionOfAtom(mview->mol, n1);
2069                                         if (n2 >= 0)
2070                                                 MoleculeToggleSelectionOfAtom(mview->mol, n2);
2071                                 } else {
2072                                         if (n2 < 0) {
2073                                                 if (!MoleculeIsAtomSelected(mview->mol, n1)) {
2074                                                         /*  If this atom is not selected, select this atom and unselect others  */
2075                                                         MoleculeSelectAtom(mview->mol, n1, 0);
2076                                                 }
2077                                         } else {
2078                                                 if (!MoleculeIsBondSelected(mview->mol, n1, n2)) {
2079                                                         /*  If this bond is not selected, select this bond and unselect others */
2080                                                         MoleculeSelectAtom(mview->mol, n1, 0);
2081                                                         MoleculeSelectAtom(mview->mol, n2, 1);
2082                                                 }
2083                                         }
2084                                 }
2085                                 if (MoleculeIsAtomSelected(mview->mol, n1))
2086                                         mview->draggingMode = kMainViewDraggingSelectedAtoms;
2087                                 else mview->draggingMode = 0;
2088                         } else mview->draggingMode = kMainViewSelectingRegion;
2089                         mview->dragEndPos[0] = mousePos[0];
2090                         mview->dragEndPos[1] = mousePos[1];
2091                         mview->dragEndPos[2] = mview->dragStartPos[2];
2092                         mview->modifierFlags = flags;
2093                         break;
2094                 }
2095                 case kTrackballCreateMode: {
2096                         if (md_is_running(mview->mol->arena)) {
2097                                 MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
2098                                 mview->draggingMode = 0;
2099                                 break;
2100                         }
2101                         /*  Draw a new bond  */
2102                         MoleculeSetSelection(mview->mol, NULL);
2103                         if (found && n2 < 0) {
2104                                 /*  An atom under mouse: create a new atom and a bond  */
2105                                 ap = ATOM_AT_INDEX(mview->mol->atoms, n1);
2106                                 mview->tempAtoms[0] = n1;
2107                                 mview->tempAtomPos[0] = ap->r;
2108                         } else {
2109                                 /*  No atom under mouse: create two new atoms and a bond  */
2110                                 mview->tempAtoms[0] = -1;
2111                                 screenPos[0] = mousePos[0];
2112                                 screenPos[1] = mousePos[1];
2113                                 screenPos[2] = mview->dragStartPos[2];
2114                                 if (MainView_convertScreenPositionToObjectPosition(mview, screenPos, objectPos) == 0) {
2115                                         mview->tempAtoms[0] = -1;  /*  Cannot create  */
2116                                         mview->draggingMode = 0;
2117                                 } else {
2118                                         mview->tempAtomPos[0].x = objectPos[0];
2119                                         mview->tempAtomPos[0].y = objectPos[1];
2120                                         mview->tempAtomPos[0].z = objectPos[2];
2121                                 /*      if (mview->mol->is_xtal_coord)
2122                                                 TransformVec(&mview->tempAtomPos[0], mview->mol->cell->rtr, &mview->tempAtomPos[0]); */
2123                                 }
2124                         }
2125                         mview->tempAtoms[1] = -1;
2126                         mview->tempAtomPos[1] = mview->tempAtomPos[0];
2127                         mview->draggingMode = kMainViewCreatingBond;
2128                         break;
2129                 }
2130                 default:
2131                         mview->draggingMode = 0;
2132                         break;
2133         }
2134         if (found)
2135                 showAtomsInInfoText(mview, n1, n2);
2136         else
2137                 showAtomsInInfoText(mview, -1, -1);
2138 /*      MainViewCallback_setNeedsDisplay(mview, 1); */
2139 }
2140
2141 void
2142 MainView_mouseDragged(MainView *mview, const float *mousePos, int flags)
2143 {
2144         float p[2];
2145         if (mview->isDragging == 0) {
2146                 if (abs(mousePos[0] - mview->dragStartPos[0]) >= 3 || abs(mousePos[1] - mview->dragStartPos[1]) >= 3)
2147                         mview->isDragging = 1;
2148                 else return;
2149         }
2150         mousePosToTrackballPos(mview, mousePos, p);
2151         switch (mview->draggingMode) {
2152                 case kMainViewMovingTrackball:
2153                         TrackballDrag(mview->track, p);
2154                         break;
2155                 case kMainViewSelectingRegion:
2156                 case kMainViewDraggingSelectedAtoms:
2157                         mview->dragEndPos[0] = mousePos[0];
2158                         mview->dragEndPos[1] = mousePos[1];
2159                         mview->dragEndPos[2] = mview->dragStartPos[2];
2160                         mview->modifierFlags = flags;
2161                         MainViewCallback_display(mview);
2162                         break;
2163                 case kMainViewCreatingBond: {
2164                         int n1, n2;
2165                         if (MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 0, 0) && n2 < 0) {
2166                                 /*  An atom under mouse  */
2167                                 Atom *ap = ATOM_AT_INDEX(mview->mol->atoms, n1);
2168                         /*      printf("n1=%d, n2=%d\n", n1, n2); */
2169                                 mview->tempAtoms[1] = n1;
2170                                 mview->tempAtomPos[1] = ap->r;
2171                         } else {
2172                                 float screenPos[3];
2173                                 double objectPos[3];
2174                                 Vector r1;
2175                                 mview->tempAtoms[1] = -1;
2176                                 /*  Convert the position of temporary atom 0 to screen position */
2177                                 r1 = mview->tempAtomPos[0];
2178                         /*      if (mview->mol->is_xtal_coord)
2179                                         TransformVec(&r1, mview->mol->cell->tr, &r1); */
2180                                 objectPos[0] = r1.x;
2181                                 objectPos[1] = r1.y;
2182                                 objectPos[2] = r1.z;
2183                                 if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, screenPos) == 0)
2184                                         break;  /*  Do nothing  */
2185                                 /*  Convert the mouse position to object position, while using the same Z depth calculated from the temporary atom 0 position  */
2186                                 screenPos[0] = mousePos[0];
2187                                 screenPos[1] = mousePos[1];
2188                                 if (MainView_convertScreenPositionToObjectPosition(mview, screenPos, objectPos) == 0)
2189                                         break;  /*  Do nothing  */
2190                                 r1.x = objectPos[0];
2191                                 r1.y = objectPos[1];
2192                                 r1.z = objectPos[2];
2193                         /*      if (mview->mol->is_xtal_coord)
2194                                         TransformVec(&r1, mview->mol->cell->rtr, &r1); */
2195                                 mview->tempAtomPos[1] = r1;
2196                         }
2197                         if (mview->tempAtoms[0] < 0)
2198                                 showAtomsInInfoText(mview, mview->tempAtoms[1], -1);
2199                         else
2200                                 showAtomsInInfoText(mview, mview->tempAtoms[0], mview->tempAtoms[1]);
2201                         MainViewCallback_display(mview);
2202                         break;
2203                 }
2204                 default:
2205                         return;
2206         }
2207         MainViewCallback_display(mview);
2208 }
2209
2210 void
2211 MainView_mouseMoved(MainView *mview, const float *mousePos, int flags)
2212 {
2213         int n1, n2, found;
2214         if (mview->isDragging)
2215                 return;
2216         found = MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 1, 0);
2217         if (found)
2218                 showAtomsInInfoText(mview, n1, n2);
2219         else
2220                 showAtomsInInfoText(mview, -1, -1);
2221 }
2222
2223 static int
2224 sCreateNewAtom(Molecule *mol, Vector pos)
2225 {
2226         Int i, j;
2227         Atom a;
2228         char name[6];
2229         memset(&a, 0, sizeof(Atom));
2230         a.occupancy = 1.0;
2231         for (i = 0; i < 1000; i++) {
2232                 snprintf(name, sizeof name, "C%03d", i);
2233                 for (j = 0; j < mol->natoms; j++) {
2234                         if (strncmp(ATOM_AT_INDEX(mol->atoms, j)->aname, name, 4) == 0)
2235                                 break;
2236                 }
2237                 if (j == mol->natoms)
2238                         break;
2239         }
2240         strncpy(a.aname, name, 4);
2241         a.atomicNumber = 6;
2242         strcpy(a.element, "C");
2243         a.type = AtomTypeEncodeToUInt("c3");
2244         a.weight = WeightForAtomicNumber(a.atomicNumber);
2245         a.r = pos;
2246         if (MolActionCreateAndPerform(mol, gMolActionAddAnAtom, &a, -1, &i) == 0)
2247                 return mol->natoms - 1;
2248         else return -1;
2249 /*      return MoleculeAddAtom(mol, &a); */
2250 }
2251
2252 void
2253 MainView_mouseUp(MainView *mview, const float *mousePos, int flags, int clickCount)
2254 {
2255         float p[2];
2256         char buf[1024];
2257
2258         mousePosToTrackballPos(mview, mousePos, p);
2259
2260         if (clickCount == 2 && mview->isDragging == 0) {
2261                 /*  Create a new molecular fragment  */
2262                 int n1, n2, found;
2263                 found = MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 0, 0);
2264                 buf[0] = 0;
2265                 if (MyAppCallback_getTextWithPrompt("Enter formula (e.g. CH2OCH3)", buf, sizeof buf) == 0)
2266                         return;
2267                 MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("si"), "cmd_fuse_or_dock", buf, n1);
2268                 goto exit;
2269         }
2270
2271         if (mview->mode == kTrackballEraseMode) {
2272                 int n1, n2, found;
2273                 found = MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 0, 0);
2274                 if (found && n1 == mview->clickedAtoms[0] && n2 == mview->clickedAtoms[1]) {
2275                         if (n2 < 0) {
2276                         /*      IntGroup *sel, *newsel;
2277                                 sel = MoleculeGetSelection(mview->mol);
2278                                 if (sel != NULL) {
2279                                         newsel = MoleculeModifySelectionByRemovingAtoms(mview->mol, sel, IntGroupNewWithPoints(n1, 1, -1));
2280                                         if (!IntGroupIsEqual(sel, newsel))
2281                                                 MolActionCreateAndPerform(mview->mol, gMolActionSetSelection, newsel);
2282                                         if (newsel != NULL)
2283                                                 IntGroupRelease(newsel);
2284                                 } */
2285                                 MolActionCreateAndPerform(mview->mol, gMolActionDeleteAnAtom, (Int)n1);
2286
2287                         } else {
2288                                 Int bn = MoleculeLookupBond(mview->mol, n1, n2);
2289                                 if (bn >= 0) {
2290                                         IntGroup *ig = IntGroupNewWithPoints(bn, 1, -1);
2291                                         MolActionCreateAndPerform(mview->mol, gMolActionDeleteBonds, ig);
2292                                         IntGroupRelease(ig);
2293                                 } else {
2294                                         fprintf(stderr, "Internal error %s:%d: bond to delete is not found", __FILE__, __LINE__);
2295                                 }
2296                         }
2297                 }
2298                 goto exit;
2299         }
2300
2301         switch (mview->draggingMode) {
2302                 case kMainViewMovingTrackball:
2303                         TrackballEndDragging(mview->track, p);
2304                         break;
2305                 case kMainViewDraggingSelectedAtoms: {
2306                         Vector offset;
2307                         if (mview->isDragging) {
2308                                 IntGroup *dupSelection = IntGroupNewFromIntGroup(mview->mol->selection);
2309                                 calcDragOffset(mview, &offset);
2310                         /*      if (mview->mol->is_xtal_coord)
2311                                         TransformVec(&offset, mview->mol->cell->rtr, &offset); */
2312                                 MolActionCreateAndPerform(mview->mol, gMolActionTranslateAtoms, &offset, dupSelection);
2313                                 IntGroupRelease(dupSelection);
2314                 /*      } else if (clickCount == 2) {
2315                                 buf[0] = 0;
2316                                 if (MyAppCallback_getTextWithPrompt("Enter formula (e.g. CH2OCH3)", buf, sizeof buf) == 0)
2317                                         return;
2318                                 MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("s"), "dock_formula", buf);
2319                 */
2320                         }
2321                         break;
2322                 }
2323                 case kMainViewSelectingRegion: {
2324                         char *selectFlags = temporarySelection(mview, mview->modifierFlags, (mview->isDragging ? 0 : 1), 0);
2325                         if (selectFlags != NULL) {
2326                                 IntGroup *ig = IntGroupNew();
2327                                 if (ig != NULL) {
2328                                         int i, natoms;
2329                                         natoms = mview->mol->natoms;
2330                                         for (i = 0; i < natoms; i++) {
2331                                                 if (selectFlags[i])
2332                                                         IntGroupAdd(ig, i, 1);
2333                                         }
2334                                         MoleculeSetSelection(mview->mol, ig);
2335                                 /*      printf("current selection = {");
2336                                         for (i = j = 0; i < natoms; i++) {
2337                                                 if (selectFlags[i]) {
2338                                                         printf("%s%.4s", (j == 0 ? "" : " "), mview->mol->atoms[i].name);
2339                                                         j++;
2340                                                 }
2341                                         }
2342                                         printf("}\n"); */
2343                                         IntGroupRelease(ig);
2344                                 }
2345                                 free(selectFlags);
2346                         }
2347                         break;
2348                 }
2349                 case kMainViewCreatingBond: {
2350                         Int b[3];
2351                         int n1, n2;
2352                         if (mview->tempAtoms[0] >= 0 && mview->tempAtoms[1] >= 0) {
2353                                 n1 = mview->tempAtoms[0];
2354                                 n2 = mview->tempAtoms[1];
2355                         } else {
2356                                 if (mview->tempAtoms[0] < 0) {
2357                                         /*  Create the first atom  */
2358                                         n1 = sCreateNewAtom(mview->mol, mview->tempAtomPos[0]);
2359                                 } else n1 = mview->tempAtoms[0];
2360                                 if (mview->tempAtoms[1] < 0) {
2361                                         /*  Create the second atom, if not too close  */
2362                                         Vector dr = ATOM_AT_INDEX(mview->mol->atoms, n1)->r;
2363                                         VecDec(dr, mview->tempAtomPos[1]);
2364                                 /*      if (mview->mol->is_xtal_coord)
2365                                                 TransformVec(&dr, mview->mol->cell->tr, &dr); */
2366                                         if (VecLength2(dr) > 0.01)
2367                                                 n2 = sCreateNewAtom(mview->mol, mview->tempAtomPos[1]);
2368                                         else n2 = -1;
2369                                 } else n2 = mview->tempAtoms[1];
2370                         }
2371                         if (n1 >= 0 && n2 >= 0 && n1 != n2) {
2372                                 b[0] = n1;
2373                                 b[1] = n2;
2374                                 b[2] = kInvalidIndex;
2375                                 MolActionCreateAndPerform(mview->mol, gMolActionAddBonds, 2, b, NULL);
2376                         /*      MoleculeAddBonds(mview->mol, b, NULL); */
2377                         }
2378                         break;
2379                 }
2380         }
2381   exit:
2382         mview->draggingMode = 0;
2383         mview->isDragging = 0;
2384         MainViewCallback_setNeedsDisplay(mview, 1);
2385         MainViewCallback_setKeyboardFocus(mview);
2386 }
2387
2388 void
2389 MainView_rotateBySlider(MainView *mview, float angle, int mode, int mouseStatus, int modifierFlags)
2390 {
2391         int i, n1, n2;
2392
2393         if (mouseStatus == 1) {  /*  mouseDown  */
2394
2395                 if (mode == kSliderRotateBondMode) {
2396                         if (MoleculeIsFragmentRotatable(mview->mol, mview->mol->selection, &n1, &n2, &(mview->rotateFragment))) {
2397                                 mview->rotateCenter = ATOM_AT_INDEX(mview->mol->atoms, n1)->r;
2398                                 mview->rotateAxis = ATOM_AT_INDEX(mview->mol->atoms, n2)->r;
2399                                 VecDec(mview->rotateAxis, mview->rotateCenter);
2400                                 if (VecLength2(mview->rotateAxis) < 1e-10)
2401                                         return;
2402                                 NormalizeVec(&(mview->rotateAxis), &(mview->rotateAxis));
2403                                 if (modifierFlags & kAltKeyMask) {
2404                                         /*  Option key: reverse the fragment  */
2405                                         IntGroupRelease(mview->rotateFragment);
2406                                         mview->rotateFragment = MoleculeFragmentExcludingAtoms(mview->mol, n2, 1, &n1);
2407                                         VecScaleSelf(mview->rotateAxis, -1);
2408                                 }
2409                         } else return;
2410                 } else if ((modifierFlags & kAltKeyMask) != 0) {
2411                         /*  Rotate selection  */
2412                         double tr[3];
2413                         if (mview->mol->selection == NULL)
2414                                 return;
2415                         mview->rotateFragment = IntGroupNewFromIntGroup(mview->mol->selection);
2416                         TrackballGetTranslate(mview->track, tr);
2417                         mview->rotateCenter.x = -mview->dimension * tr[0];
2418                         mview->rotateCenter.y = -mview->dimension * tr[1];
2419                         mview->rotateCenter.z = -mview->dimension * tr[2];
2420                         if (mode == kSliderRotateXMode) {
2421                                 VecCross(mview->rotateAxis, mview->up, mview->lookto);
2422                         } else {
2423                                 mview->rotateAxis = mview->up;
2424                                 VecScaleSelf(mview->rotateAxis, -1);
2425                         }
2426                 } else {
2427                         /*  Rotate along x/y axis (no coordinate transform)  */
2428                         TrackballStartDragging(mview->track, NULL, kTrackballRotateMode);
2429                         mview->rotateFragment = NULL;  /*  This is probably not necessary  */
2430                 }
2431                 if (mview->rotateFragment != NULL) {
2432                         /*  Save the original position  */
2433                         n1 = IntGroupGetCount(mview->rotateFragment);
2434                         mview->rotateFragmentOldPos = (Vector *)calloc(sizeof(Vector), n1);
2435                         if (mview->rotateFragmentOldPos == NULL) {
2436                                 IntGroupRelease(mview->rotateFragment);
2437                                 mview->rotateFragment = NULL;
2438                                 return;
2439                         }
2440                         for (i = 0; i < n1; i++) {
2441                                 n2 = IntGroupGetNthPoint(mview->rotateFragment, i);
2442                                 mview->rotateFragmentOldPos[i] = (ATOM_AT_INDEX(mview->mol->atoms, n2))->r;
2443                         }
2444                 }
2445
2446         } else {  /*  mouseDragged or mouseUp */
2447
2448                 if (mview->rotateFragment != NULL) {
2449
2450                         /*  Restore original positions  */
2451                         n1 = IntGroupGetCount(mview->rotateFragment);
2452                         for (i = 0; i < n1; i++) {
2453                                 n2 = IntGroupGetNthPoint(mview->rotateFragment, i);
2454                                 (ATOM_AT_INDEX(mview->mol->atoms, n2))->r = mview->rotateFragmentOldPos[i];
2455                         }
2456                         /*  Rotate  */
2457                         if (mouseStatus == 2) { /* dragged */
2458                                 MoleculeRotate(mview->mol, &(mview->rotateAxis), angle, &(mview->rotateCenter), mview->rotateFragment);
2459                         } else {  /* mouse up */
2460                                 IntGroup *ig = IntGroupNewFromIntGroup(mview->rotateFragment);
2461                                 MolActionCreateAndPerform(mview->mol, gMolActionRotateAtoms, &(mview->rotateAxis), (double)angle, &(mview->rotateCenter), ig);
2462                                 IntGroupRelease(ig);
2463                                 IntGroupRelease(mview->rotateFragment);
2464                                 mview->rotateFragment = NULL;
2465                                 free(mview->rotateFragmentOldPos);
2466                                 mview->rotateFragmentOldPos = NULL;
2467                         }
2468
2469                 } else {
2470                         double quat[4];
2471                         double cs = cos(angle / 2);
2472                         double sn = sin(angle / 2);
2473                         if (mouseStatus == 0) { /* mouseUp */
2474                                 TrackballEndDragging(mview->track, NULL);
2475                         } else { /* mouseDragged */
2476                                 quat[0] = cs;
2477                                 if (mode == kSliderRotateXMode) {
2478                                         /*  Rotation along x axis  */
2479                                         quat[1] = -sn;
2480                                         quat[2] = quat[3] = 0;
2481                                 } else {
2482                                         quat[2] = sn;
2483                                         quat[1] = quat[3] = 0;
2484                                 }
2485                                 TrackballSetTemporaryRotation(mview->track, quat);
2486                         }
2487                 }
2488         }
2489
2490         MainViewCallback_setNeedsDisplay(mview, 1);
2491         MainViewCallback_setKeyboardFocus(mview);
2492 }
2493
2494 #pragma mark ====== Menu Commands ======
2495
2496 void
2497 MainView_selectAll(MainView *mview)
2498 {
2499         IntGroup *ig;
2500         if (mview == NULL || mview->mol == NULL)
2501                 return;
2502         ig = IntGroupNew();
2503         if (ig == NULL)
2504                 return;
2505         IntGroupAdd(ig, 0, mview->mol->natoms);
2506         MoleculeSetSelection(mview->mol, ig);
2507         IntGroupRelease(ig);
2508         MainViewCallback_setNeedsDisplay(mview, 1);
2509 }
2510
2511 void
2512 MainView_selectFragment(MainView *mview)
2513 {
2514         IntGroup *ig;
2515         if (mview == NULL || mview->mol == NULL)
2516                 return;
2517         ig = MoleculeFragmentWithAtomGroups(mview->mol, MoleculeGetSelection(mview->mol), NULL);
2518         if (ig == NULL)
2519                 return;
2520         MoleculeSetSelection(mview->mol, ig);
2521         IntGroupRelease(ig);
2522         MainViewCallback_setNeedsDisplay(mview, 1);
2523 }
2524
2525 void
2526 MainView_selectReverse(MainView *mview)
2527 {
2528         IntGroup *ig;
2529         if (mview == NULL || mview->mol == NULL)
2530                 return;
2531         ig = IntGroupNewFromIntGroup(MoleculeGetSelection(mview->mol));
2532         IntGroupReverse(ig, 0, mview->mol->natoms);
2533         MoleculeSetSelection(mview->mol, ig);
2534         IntGroupRelease(ig);
2535         MainViewCallback_setNeedsDisplay(mview, 1);
2536 }
2537
2538 void
2539 MainView_centerSelection(MainView *mview)
2540 {
2541         IntGroup *ig;
2542         Vector c;
2543         double tr[3];
2544         if (mview == NULL || mview->mol == NULL)
2545                 return;
2546         ig = MoleculeGetSelection(mview->mol);
2547         MoleculeCenterOfMass(mview->mol, &c, ig);
2548         tr[0] = -c.x / mview->dimension;
2549         tr[1] = -c.y / mview->dimension;
2550         tr[2] = -c.z / mview->dimension;
2551         TrackballSetTranslate(mview->track, tr);
2552         MainViewCallback_setNeedsDisplay(mview, 1);
2553 }
2554
2555 #pragma mark ====== Pasteboard Support ======
2556
2557 int
2558 MainView_copy(MainView *mview)
2559 {
2560         Molecule *mol, *mol2;
2561         IntGroup *sel;
2562         Int len, result, time;
2563         char *p;
2564         if (mview == NULL || (mol = mview->mol) == NULL || (sel = MoleculeGetSelection(mol)) == NULL)
2565                 return 1;
2566         if (MoleculeExtract(mol, &mol2, MoleculeGetSelection(mol), 1) != 0
2567         || mol2 == NULL
2568         || (p = MoleculeSerialize(mol2, &len, &time)) == NULL)
2569                 return 2;
2570         result = MoleculeCallback_writeToPasteboard(gMoleculePasteboardType, p, len);
2571         if (result != 0)
2572                 return result;
2573         if (mol2 != NULL)
2574                 MoleculeRelease(mol2);
2575         mview->pasteTimeStamp = time;
2576         mview->pasteCount = 0;
2577         return 0;
2578 }
2579
2580 int
2581 MainView_delete(MainView *mview)
2582 {
2583         int result;
2584         Molecule *mol = mview->mol;
2585         if (md_is_running(mview->mol->arena)) {
2586                 MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
2587                 return -1;
2588         }
2589         result = MolActionCreateAndPerform(mol, gMolActionUnmergeMolecule, MoleculeGetSelection(mol));
2590         if (result == 0)
2591                 result = MolActionCreateAndPerform(mol, gMolActionSetSelection, NULL);
2592         return result;
2593 }
2594
2595 int
2596 MainView_cut(MainView *mview)
2597 {
2598         int result;
2599         if (md_is_running(mview->mol->arena)) {
2600                 MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
2601                 return -1;
2602         }
2603         result = MainView_copy(mview);
2604         if (result == 0)
2605                 result = MainView_delete(mview);
2606         return result;
2607 }
2608
2609 int
2610 MainView_isPastable(MainView *mview)
2611 {
2612         if (MoleculeCallback_isDataInPasteboard(gMoleculePasteboardType))
2613                 return 1;
2614         else return 0;
2615 }
2616
2617 int
2618 MainView_paste(MainView *mview)
2619 {
2620         void *p;
2621         Int len, result, time;
2622         Molecule *mol2;
2623         IntGroup *sel;
2624
2625         if (md_is_running(mview->mol->arena)) {
2626                 MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
2627                 return -1;
2628         }
2629         if (!MainView_isPastable(mview))
2630                 return 1;
2631         if (MoleculeCallback_readFromPasteboard(gMoleculePasteboardType, &p, &len) != 0)
2632                 return 2;
2633         if ((mol2 = MoleculeDeserialize(p, len, &time)) == NULL) {
2634                 free(p);
2635                 return 3;
2636         }
2637         free(p);
2638
2639         if (time == mview->pasteTimeStamp) {
2640                 /*  Offset the pasted fragment by something  */
2641                 Vector v;
2642                 mview->pasteCount++;
2643                 v.x = v.y = v.z = mview->pasteCount * 0.5 + sin((double)mview->pasteCount) * 0.1; /* sin(..) is to avoid accidental coincidence  */
2644                 MoleculeTranslate(mol2, &v, NULL);
2645         } else {
2646                 mview->pasteTimeStamp = time;
2647                 mview->pasteCount = 0;
2648         }
2649
2650         sel = MoleculeGetSelection(mview->mol);
2651         if (sel == NULL || IntGroupGetCount(sel) == 0) {
2652                 /*  Remove dummy atoms from mol2  */
2653                 Atom *ap;
2654                 int i;
2655                 sel = IntGroupNew();
2656                 for (i = 0, ap = mol2->atoms; i < mol2->natoms; i++, ap = ATOM_NEXT(ap)) {
2657                         if (ap->atomicNumber == 0 && ap->aname[0] == '_') {
2658                                 /*  Dummy atom  */
2659                                 IntGroupAdd(sel, i, 1);
2660                         }
2661                 }
2662                 if (IntGroupGetCount(sel) > 0)
2663                         MoleculeUnmerge(mol2, NULL, sel, 0, NULL, NULL, 0);
2664                 IntGroupRelease(sel);
2665         }
2666
2667         result = MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("M"), "dock_fragment", mol2);
2668         MoleculeRelease(mol2);
2669         return result;
2670 }
2671
2672 static int
2673 sMainView_AppendParameterToSerialBuffer(void **pb_ptr, Int *pb_len, Int parType, Int count, const UnionPar *up)
2674 {
2675         int len = sizeof(Int) * 2 + sizeof(UnionPar) * count;
2676         char *p;
2677         if (*pb_len == 0) {
2678                 *pb_ptr = malloc(len + 1);
2679                 if (*pb_ptr == NULL)
2680                         return -1;
2681                 p = (char *)(*pb_ptr);
2682                 *pb_len = len + 1;
2683         } else {
2684                 *pb_ptr = realloc(*pb_ptr, *pb_len + len);
2685                 if (*pb_ptr == NULL)
2686                         return -1;
2687                 p = (char *)(*pb_ptr) + *pb_len - 1;  /*  *pb_len includes the end mark  */
2688                 *pb_len += len;
2689         }
2690         *((Int *)p) = 0;
2691         *p = parType;
2692         p += sizeof(Int);
2693         *((Int *)p) = count;
2694         p += sizeof(Int);
2695         memmove(p, up, sizeof(UnionPar) * count);
2696         p += sizeof(UnionPar) * count;
2697         *p = 0;  /*  End mark  */
2698         return 0;
2699 }
2700
2701 int
2702 MainView_pasteParameters(MainView *mview)
2703 {
2704         char *p, *pp;
2705         Int len, i;
2706         IntGroup *newsel;
2707         if (mview == NULL || mview->mol == NULL)
2708                 return -1;
2709         if (mview->mol->par == NULL)
2710                 mview->mol->par = ParameterNew();
2711         if (!MoleculeCallback_isDataInPasteboard(gParameterPasteboardType))
2712                 return 1;
2713         if (MoleculeCallback_readFromPasteboard(gParameterPasteboardType, (void **)&p, &len) != 0)
2714                 return 2;
2715         pp = p;
2716         newsel = IntGroupNew();
2717         while (*p != 0) {
2718                 int count, c;
2719                 UnionPar *up;
2720                 IntGroup *ig;
2721                 int parType = *p;
2722                 if (parType < kFirstParType || parType > kLastParType)
2723                         break;
2724                 p += sizeof(Int);
2725                 count = *((Int *)p);
2726                 p += sizeof(Int);
2727                 up = (UnionPar *)calloc(sizeof(UnionPar), count);
2728                 memmove(up, p, sizeof(UnionPar) * count);
2729
2730                 /*  The global parameters become local when pasted  */
2731                 for (i = 0; i < count; i++) {
2732                         if (up[i].bond.src > 0)
2733                                 up[i].bond.src = 0;
2734                 }
2735
2736                 c = ParameterGetCountForType(mview->mol->par, parType);
2737                 ig = IntGroupNewWithPoints(c, count, -1);
2738                 MolActionCreateAndPerform(mview->mol, gMolActionAddParameters, parType, ig, count, up);
2739                 free(up);
2740                 IntGroupRelease(ig);
2741                 p += sizeof(UnionPar) * count;
2742                 c = ParameterTableGetRowFromTypeAndIndex(mview->mol->par, parType, c);
2743                 IntGroupAdd(newsel, c, count);
2744         }
2745         free(pp);
2746
2747         /*  Select newly pasted parameters  */
2748         MainViewCallback_setTableSelection(mview, newsel);
2749         IntGroupRelease(newsel);
2750
2751         /*  Request MD rebuild  */
2752         mview->mol->needsMDRebuild = 1;
2753
2754         /*  Suppress clear of parameter table selection  */
2755         mview->mol->parameterTableSelectionNeedsClear = 0;
2756
2757         MoleculeCallback_notifyModification(mview->mol, 0);
2758         return 0;
2759 }
2760
2761 /*  flags: 1, delete; 2, copy; 3, cut  */
2762 int
2763 MainView_copyOrCutParameters(MainView *mview, int flags)
2764 {
2765         IntGroup *ig = MainViewCallback_getTableSelection(mview);
2766         IntGroup *ig2;
2767         int i, n, idx, type = -1, t1;
2768         Parameter *par = mview->mol->par;
2769         void *pb_ptr = NULL;
2770         Int pb_len = 0;
2771         if (ig == NULL || (i = IntGroupGetCount(ig)) == 0)
2772                 return 0;
2773         i--;
2774         ig2 = NULL;
2775         while (1) {
2776                 n = IntGroupGetNthPoint(ig, i);
2777                 if (n >= 0)
2778                         idx = ParameterTableGetItemIndex(par, n, &t1);
2779                 else {
2780                         idx = -1;
2781                         t1 = -2;
2782                 }
2783                 if (t1 != type) {
2784                         /*  Process Parameters for the last group  */
2785                         if (type >= kFirstParType && ig2 != NULL && (n = IntGroupGetCount(ig2)) > 0) {
2786                                 UnionPar *up = (UnionPar *)calloc(sizeof(UnionPar), n);
2787                                 if (flags & 1) {
2788                                         MolAction *act;
2789                                         if (ParameterDelete(par, type, up, ig2) < 0)
2790                                                 return -1;
2791                                         act = MolActionNew(gMolActionAddParameters, type, ig2, n, up);
2792                                         MolActionCallback_registerUndo(mview->mol, act);
2793                                         MolActionRelease(act);
2794                                 } else {
2795                                         if (ParameterCopy(par, type, up, ig2) < 0)
2796                                                 return -1;
2797                                 }
2798                                 if (flags & 2) {
2799                                         if (sMainView_AppendParameterToSerialBuffer(&pb_ptr, &pb_len, type, n, up) != 0)
2800                                                 return -1;
2801                                 }
2802                                 free(up);
2803                                 IntGroupRelease(ig2);
2804                                 ig2 = NULL;
2805                         }
2806                         if (t1 == -2)
2807                                 break;
2808                         type = t1;
2809                 }
2810                 if (idx >= 0) {
2811                         if (ig2 == NULL)
2812                                 ig2 = IntGroupNew();
2813                         IntGroupAdd(ig2, idx, 1);
2814                 }
2815                 i--;
2816         }
2817         if (ig2 != NULL)
2818                 IntGroupRelease(ig2);
2819         IntGroupRelease(ig);
2820
2821         if (flags & 2) {
2822                 n = MoleculeCallback_writeToPasteboard(gParameterPasteboardType, pb_ptr, pb_len);
2823                 if (n != 0)
2824                         return n;
2825         }
2826         if (flags & 1) {
2827                 /*  Clear selection  */
2828                 MainViewCallback_setTableSelection(mview, NULL);
2829                 /*  Request MD rebuild  */
2830                 mview->mol->needsMDRebuild = 1;
2831         }
2832         MoleculeCallback_notifyModification(mview->mol, 0);
2833         return 0;
2834 }
2835
2836 #pragma mark ====== Table View (also for gBuiltinParameters) ======
2837
2838 /*  As a special case, the table view functions will handle gBuiltinParameters when mview == NULL.  */
2839
2840 typedef struct ColumnInfoRecord {
2841         char *name;
2842         int width;
2843         int editable;
2844 } ColumnInfoRecord;
2845
2846 static ColumnInfoRecord sAtomColumns[] = {
2847         {"atom", 4, 0}, {"name", 4, 1}, {"type", 4, 1}, {"element", 4, 1}, {"residue", 6, 1},
2848         {"x", 6, 1}, {"y", 6, 1}, {"z", 6, 1}, {"charge", 6, 1}, {NULL}
2849 };
2850 static ColumnInfoRecord sBondColumns[] = {
2851         {"atoms", 9, 0}, {"names", 9, 0}, {"type", 9, 0}, {"length", 8, 0},
2852         {"par type", 8, 0}, {"k", 8, 0}, {"r0", 8, 0}, {NULL}
2853 };
2854 static ColumnInfoRecord sAngleColumns[] = {
2855         {"atoms", 12, 0}, {"names", 12, 0}, {"type", 12, 0}, {"angle", 8, 0},
2856         {"par type", 8, 0}, {"k", 8, 0}, {"a0", 8, 0}, {NULL}
2857 };
2858 static ColumnInfoRecord sDihedralColumns[] = {
2859         {"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"dihedral", 8, 0},
2860         {"par type", 8, 0}, {"k", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
2861 };
2862 static ColumnInfoRecord sImproperColumns[] = {
2863         {"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"improper", 8, 0},
2864         {"par type", 8, 0}, {"k", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
2865 };
2866 static ColumnInfoRecord sParameterColumns[] = {
2867         {"class", 5, 0}, {"type", 9, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0},
2868         {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"src", 8, 0}, {"comment", 25, 0}, {NULL}
2869 };
2870 static ColumnInfoRecord sUnitCellColumns[] = {
2871         {"name", 6, 0}, {"values", 9, 1}, {"", 9, 1}, {"", 9, 1}, {NULL}
2872 };
2873 static ColumnInfoRecord sXtalCoordColumns[] = {
2874         {"atom", 4, 0}, {"name", 4, 1}, {"element", 4, 1},
2875         {"frx", 6, 1}, {"fry", 6, 1}, {"frz", 6, 1},
2876         {"occ", 6, 1}, {"temp", 6, 1}, {NULL}
2877 };
2878 static ColumnInfoRecord sMOEnergyColumns[] = {
2879         {"MO", 5, 0}, {"alpha energy", 12, 0}, {"beta energy", 12, 0}, {NULL}
2880 };
2881 static ColumnInfoRecord *sColumnInfo[] = {
2882         sAtomColumns, sBondColumns, sAngleColumns, sDihedralColumns,
2883         sImproperColumns, NULL, sParameterColumns, NULL,
2884         sUnitCellColumns, sXtalCoordColumns, NULL, sMOEnergyColumns
2885 };
2886 static char *sTableTitles[] = {
2887         "atoms", "bonds", "angles", "dihedrals", "impropers", "--------",
2888         "MM/MD pars", "--------", "unit cell", "xtal coords", "--------", "MO energy"
2889 };
2890
2891 void
2892 MainView_tableTitleForIndex(MainView *mview, int idx, char *buf, int bufsize)
2893 {
2894         if (mview == NULL)
2895                 idx = kMainViewParameterTableIndex;
2896         if (idx < 0 || idx >= sizeof(sColumnInfo) / sizeof(sColumnInfo[0])) {
2897                 buf[0] = 0;
2898                 return;
2899         }
2900         snprintf(buf, bufsize, "%s", sTableTitles[idx]);
2901 }
2902
2903 int
2904 MainView_createColumnsForTableAtIndex(MainView *mview, int idx)
2905 {
2906         int i;
2907
2908         if (mview == NULL)
2909                 idx = kMainViewParameterTableIndex;
2910         if (idx < 0 || idx >= sizeof(sColumnInfo) / sizeof(sColumnInfo[0]))
2911                 return 0;  /*  Invalid index  */
2912         if (sColumnInfo[idx] == NULL)
2913                 return 0;  /*  Invalid index  */
2914
2915         /*  Remove all existing columns  */
2916         while (MainViewCallback_removeTableColumnAtIndex(mview, 0) > 0);
2917
2918         if (mview != NULL)
2919                 mview->tableIndex = idx;
2920
2921         /*  Create columns  */
2922         for (i = 0; ; i++) {
2923                 int width;
2924                 ColumnInfoRecord *recp = &(sColumnInfo[idx][i]);
2925                 if (recp->name == NULL)
2926                         break;
2927                 width = recp->width;
2928                 if (mview == 0 && strcmp(recp->name, "comment") == 0)
2929                         width = 80;
2930                 MainViewCallback_addTableColumn(mview, recp->name, width, recp->editable);
2931         }
2932
2933         return 1;
2934 }
2935
2936 void
2937 MainView_refreshTable(MainView *mview)
2938 {
2939         /*  Reload data  */
2940         if (mview != NULL && mview->mol != NULL) {
2941                 MainView_refreshCachedInfo(mview);
2942                 if (mview->tableIndex == kMainViewBondTableIndex ||
2943                         mview->tableIndex == kMainViewAngleTableIndex ||
2944                         mview->tableIndex == kMainViewDihedralTableIndex ||
2945                         mview->tableIndex == kMainViewImproperTableIndex ||
2946                         mview->tableIndex == kMainViewParameterTableIndex) {
2947                         /*  Check the parameter table  */
2948                         if (mview->mol->arena == NULL)
2949                                 md_arena_new(mview->mol);
2950                         if (mview->mol->needsMDRebuild || mview->mol->arena->par == NULL) {
2951                                 /*  Here we do not call MoleculePrepareMDArena(), because we do not want
2952                                         to modify Molecule by just redrawing tables.
2953                                     (But MoleculePrepareMDArena() *is* called when the parameter
2954                                         table is selected.  */
2955                                 md_prepare(mview->mol->arena, 1);
2956                         }
2957                 }
2958         }
2959
2960         MainViewCallback_reloadTableData(mview);
2961
2962         if (mview != NULL && mview->mol != NULL && mview->tableIndex >= kMainViewAtomTableIndex && mview->tableIndex <= kMainViewImproperTableIndex)
2963                 MainViewCallback_setTableSelection(mview, mview->tableSelection);
2964 }
2965
2966 int
2967 MainView_numberOfRowsInTable(MainView *mview)
2968 {
2969         if (mview == NULL)
2970                 return ParameterTableNumberOfRows(gBuiltinParameters);
2971         if (mview->mol == NULL)
2972                 return 0;
2973         switch (mview->tableIndex) {
2974                 case kMainViewParameterTableIndex:
2975                         return ParameterTableNumberOfRows(mview->mol->par);
2976                 case kMainViewUnitCellTableIndex:
2977                         return 13; /* a, b, c, alpha, beta, gamma, a_valid, b_valid, c_valid, av, bv, cv, ov */
2978                 case kMainViewAtomTableIndex:
2979                 case kMainViewBondTableIndex:
2980                 case kMainViewAngleTableIndex:
2981                 case kMainViewDihedralTableIndex:
2982                 case kMainViewImproperTableIndex:
2983                 case kMainViewXtalCoordTableIndex:
2984                         if (mview->tableCache == NULL)
2985                                 MainView_refreshCachedInfo(mview);
2986                         return IntGroupGetCount(mview->tableCache);
2987                 case kMainViewMOTableIndex:
2988                         return (mview->mol->bset != NULL ? mview->mol->bset->ncomps : 0);
2989                 default:
2990                         return 0;
2991         }
2992 }
2993
2994 int
2995 MainView_indexToTableRow(MainView *mview, int idx)
2996 {
2997         if (mview == NULL)
2998                 return -1;
2999         switch (mview->tableIndex) {
3000                 case kMainViewAtomTableIndex:
3001                 case kMainViewBondTableIndex:
3002                 case kMainViewAngleTableIndex:
3003                 case kMainViewDihedralTableIndex:
3004                 case kMainViewImproperTableIndex:
3005                 case kMainViewXtalCoordTableIndex:
3006                         return IntGroupLookupPoint(mview->tableCache, idx);
3007                 default:
3008                         return idx;
3009         }
3010 }
3011
3012 int
3013 MainView_tableRowToIndex(MainView *mview, int row)
3014 {
3015         if (mview == NULL)
3016                 return -1;
3017         switch (mview->tableIndex) {
3018                 case kMainViewAtomTableIndex:
3019                 case kMainViewBondTableIndex:
3020                 case kMainViewAngleTableIndex:
3021                 case kMainViewDihedralTableIndex:
3022                 case kMainViewImproperTableIndex:
3023                 case kMainViewXtalCoordTableIndex:
3024                         return IntGroupGetNthPoint(mview->tableCache, row);
3025                 default:
3026                         return row;
3027         }
3028 }
3029
3030 static char *
3031 sAtomDescription(Atom *ap, char *buf, int bufsize)
3032 {
3033         snprintf(buf, bufsize, "%d:%.4s", ap->resSeq, ap->aname);
3034         return buf;
3035 }
3036
3037 static UnionPar *
3038 sParameterOfTypeAtIndex(Molecule *mol, int parType, int idx)
3039 {
3040         int i;
3041         if (mol->arena == NULL || mol->needsMDRebuild || mol->arena->par == NULL)
3042                 return NULL;
3043         switch (parType) {
3044                 case kVdwParType:
3045                         i = mol->arena->vdw_par_i[idx];
3046                         if (i >= 0 && i < mol->arena->par->nvdwPars)
3047                                 return (UnionPar *)(mol->arena->par->vdwPars + i);
3048                         else return NULL;
3049                 case kBondParType:
3050                         i = mol->arena->bond_par_i[idx];
3051                         if (i >= 0 && i < mol->arena->par->nbondPars)
3052                                 return (UnionPar *)(mol->arena->par->bondPars + i);
3053                         else return NULL;
3054                 case kAngleParType:
3055                         i = mol->arena->angle_par_i[idx];
3056                         if (i >= 0 && i < mol->arena->par->nanglePars)
3057                                 return (UnionPar *)(mol->arena->par->anglePars + i);
3058                         else return NULL;
3059                 case kDihedralParType:
3060                         i = mol->arena->dihedral_par_i[idx];
3061                         if (i >= 0 && i < mol->arena->par->ndihedralPars)
3062                                 return (UnionPar *)(mol->arena->par->dihedralPars + i);
3063                         else return NULL;
3064                 case kImproperParType:
3065                         i = mol->arena->improper_par_i[idx];
3066                         if (i >= 0 && i < mol->arena->par->nimproperPars)
3067                                 return (UnionPar *)(mol->arena->par->improperPars + i);
3068                         else return NULL;
3069         }
3070         return NULL;
3071 }
3072
3073 void
3074 MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufsize)
3075 {
3076         int idx;
3077         Int *ip;
3078         Atom *ap[4];
3079         char descbuf[4][20], typebuf[4][8];
3080         Molecule *mol;
3081
3082         if (mview == NULL)
3083                 return ParameterTableGetItemText(gBuiltinParameters, column, row, buf, bufsize);
3084
3085         if (mview == NULL || (mol = mview->mol) == NULL) {
3086                 buf[0] = 0;
3087                 return;
3088         }
3089
3090         if (mview->tableIndex == kMainViewParameterTableIndex)
3091                 return ParameterTableGetItemText(mview->mol->par, column, row, buf, bufsize);
3092
3093         if (mview->tableIndex == kMainViewAtomTableIndex) { /* Atoms */
3094                 idx = MainView_tableRowToIndex(mview, row);
3095                 ap[0] = ATOM_AT_INDEX(mol->atoms, idx);
3096                 switch (column) {
3097                         case 0: snprintf(buf, bufsize, "%d", idx); break;
3098                         case 1: snprintf(buf, bufsize, "%.4s", ap[0]->aname); break;
3099                         case 2: snprintf(buf, bufsize, "%.6s", AtomTypeDecodeToString(ap[0]->type, NULL)); break;
3100                         case 3: snprintf(buf, bufsize, "%.2s", ap[0]->element); break;
3101                         case 4:
3102                                 if (ap[0]->resSeq == 0)
3103                                         buf[0] = 0;
3104                                 else
3105                                         snprintf(buf, bufsize, "%.4s.%d", ap[0]->resName, ap[0]->resSeq);
3106                                 break;
3107                         case 5: snprintf(buf, bufsize, "%.3f", ap[0]->r.x); break;
3108                         case 6: snprintf(buf, bufsize, "%.3f", ap[0]->r.y); break;
3109                         case 7: snprintf(buf, bufsize, "%.3f", ap[0]->r.z); break;
3110                         case 8: snprintf(buf, bufsize, "%.3f", ap[0]->charge); break;
3111                         default: buf[0] = 0; break;
3112                 }
3113         } else if (mview->tableIndex == kMainViewBondTableIndex) { /* Bonds */
3114                 idx = MainView_tableRowToIndex(mview, row);
3115                 ip = mol->bonds + idx * 2;
3116                 ap[0] = ATOM_AT_INDEX(mol->atoms, ip[0]);
3117                 ap[1] = ATOM_AT_INDEX(mol->atoms, ip[1]);
3118                 switch (column) {
3119                         case 0: snprintf(buf, bufsize, "%d-%d", ip[0], ip[1]); break;
3120                         case 1: snprintf(buf, bufsize, "%s-%s", sAtomDescription(ap[0], descbuf[0], 20), sAtomDescription(ap[1], descbuf[1], 20)); break;
3121                         case 2: snprintf(buf, bufsize, "%.6s-%.6s", AtomTypeDecodeToString(ap[0]->type, typebuf[0]), AtomTypeDecodeToString(ap[1]->type, typebuf[1])); break;
3122                         case 3:
3123                                 snprintf(buf, bufsize, "%.3f", MoleculeMeasureBond(mview->mol, &(ap[0]->r), &(ap[1]->r)));
3124                                 break;
3125                         case 4:
3126                         case 5:
3127                         case 6: {
3128                                 BondPar *bp = (BondPar *)sParameterOfTypeAtIndex(mol, kBondParType, idx);
3129                                 if (bp == NULL)
3130                                         buf[0] = 0;
3131                                 else {
3132                                         if (column == 4)
3133                                                 snprintf(buf, bufsize, "%.6s-%.6s", AtomTypeDecodeToString(bp->type1, typebuf[0]), AtomTypeDecodeToString(bp->type2, typebuf[1]));
3134                                         else
3135                                                 snprintf(buf, bufsize, "%.3f", (column == 5 ? bp->k * INTERNAL2KCAL : bp->r0));
3136                                 }
3137                                 break;
3138                         }
3139                         default: buf[0] = 0; break;
3140                 }
3141         } else if (mview->tableIndex == kMainViewAngleTableIndex) { /* Angles */
3142                 idx = MainView_tableRowToIndex(mview, row);
3143                 ip = mol->angles + idx * 3;
3144                 ap[0] = ATOM_AT_INDEX(mol->atoms, ip[0]);
3145                 ap[1] = ATOM_AT_INDEX(mol->atoms, ip[1]);
3146                 ap[2] = ATOM_AT_INDEX(mol->atoms, ip[2]);
3147                 switch (column) {
3148                         case 0: snprintf(buf, bufsize, "%d-%d-%d", ip[0], ip[1], ip[2]); break;
3149                         case 1: snprintf(buf, bufsize, "%s-%s-%s", sAtomDescription(ap[0], descbuf[0], 20), sAtomDescription(ap[1], descbuf[1], 20), sAtomDescription(ap[2], descbuf[2], 20)); break;
3150                         case 2: snprintf(buf, bufsize, "%.6s-%.6s-%.6s", AtomTypeDecodeToString(ap[0]->type, typebuf[0]), AtomTypeDecodeToString(ap[1]->type, typebuf[1]), AtomTypeDecodeToString(ap[2]->type, typebuf[2])); break;
3151                         case 3:
3152                                 snprintf(buf, bufsize, "%.3f", MoleculeMeasureAngle(mview->mol, &(ap[0]->r), &(ap[1]->r), &(ap[2]->r)));
3153                                 break;
3154                         case 4:
3155                         case 5:
3156                         case 6: {
3157                                 AnglePar *anp = (AnglePar *)sParameterOfTypeAtIndex(mol, kAngleParType, idx);
3158                                 if (anp == NULL)
3159                                         buf[0] = 0;
3160                                 else {
3161                                         if (column == 4)
3162                                                 snprintf(buf, bufsize, "%.6s-%.6s-%.6s", AtomTypeDecodeToString(anp->type1, typebuf[0]), AtomTypeDecodeToString(anp->type2, typebuf[1]), AtomTypeDecodeToString(anp->type3, typebuf[2]));
3163                                         else
3164                                                 snprintf(buf, bufsize, "%.3f", (column == 5 ? anp->k * INTERNAL2KCAL : anp->a0 * kRad2Deg));
3165                                 }
3166                                 break;
3167                         }
3168                         default: buf[0] = 0; break;
3169                 }
3170         } else if (mview->tableIndex == kMainViewDihedralTableIndex || mview->tableIndex == kMainViewImproperTableIndex) { /* Dihedrals, Impropers */
3171                 int f = (mview->tableIndex == kMainViewDihedralTableIndex);
3172                 idx = MainView_tableRowToIndex(mview, row);
3173                 ip = (f ? mview->mol->dihedrals : mview->mol->impropers) + idx * 4;
3174                 ap[0] = ATOM_AT_INDEX(mview->mol->atoms, ip[0]);
3175                 ap[1] = ATOM_AT_INDEX(mview->mol->atoms, ip[1]);
3176                 ap[2] = ATOM_AT_INDEX(mview->mol->atoms, ip[2]);
3177                 ap[3] = ATOM_AT_INDEX(mview->mol->atoms, ip[3]);
3178                 switch (column) {
3179                         case 0: snprintf(buf, bufsize, "%d-%d-%d-%d", ip[0], ip[1], ip[2], ip[3]); break;
3180                         case 1: snprintf(buf, bufsize, "%s-%s-%s-%s", sAtomDescription(ap[0], descbuf[0], 20), sAtomDescription(ap[1], descbuf[1], 20), sAtomDescription(ap[2], descbuf[2], 20), sAtomDescription(ap[3], descbuf[3], 20)); break;
3181                         case 2: snprintf(buf, bufsize, "%.6s-%.6s-%.6s-%.6s", AtomTypeDecodeToString(ap[0]->type, typebuf[0]), AtomTypeDecodeToString(ap[1]->type, typebuf[1]), AtomTypeDecodeToString(ap[2]->type, typebuf[2]), AtomTypeDecodeToString(ap[3]->type, typebuf[3])); break;
3182                         case 3:
3183                                 snprintf(buf, bufsize, "%.3f", MoleculeMeasureDihedral(mview->mol, &(ap[0]->r), &(ap[1]->r), &(ap[2]->r), &(ap[3]->r)));
3184                                 break;
3185                         case 4:
3186                         case 5:
3187                         case 6:
3188                         case 7: {
3189                                 TorsionPar *tp = (TorsionPar *)sParameterOfTypeAtIndex(mol, (f ? kDihedralParType : kImproperParType), idx);
3190                                 if (tp == NULL)
3191                                         buf[0] = 0;
3192                                 else if (column == 4)
3193                                         snprintf(buf, bufsize, "%.6s-%.6s-%.6s-%.6s", AtomTypeDecodeToString(tp->type1, typebuf[0]), AtomTypeDecodeToString(tp->type2, typebuf[1]), AtomTypeDecodeToString(tp->type3, typebuf[2]), AtomTypeDecodeToString(tp->type4, typebuf[3]));
3194                                 else if (column == 6)
3195                                         snprintf(buf, bufsize, "%d", tp->period[0]);
3196                                 else snprintf(buf, bufsize, "%.3f", (column == 5 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
3197                                 break;
3198                         }
3199                         default: buf[0] = 0; break;
3200                 }
3201         } else if (mview->tableIndex == kMainViewUnitCellTableIndex) { /* Unit cell info */
3202                 buf[0] = 0;
3203                 if (mview->mol->cell != NULL) {
3204                         static const char *unitCellRowTitles[] = {"a", "b", "c", "alpha", "beta", "gamma", "a_valid", "b_valid", "c_valid", "a_vec", "b_vec", "c_vec", "origin"};
3205                         if (column == 0)
3206                                 snprintf(buf, bufsize, "%s", unitCellRowTitles[row]);
3207                         else {
3208                                 switch(row) {
3209                                         case 0: case 1: case 2:
3210                                                 if (column == 1)
3211                                                         snprintf(buf, bufsize, "%.5f", mview->mol->cell->cell[row]);
3212                                                 else if (column == 2 && mview->mol->cell->has_sigma)
3213                                                         snprintf(buf, bufsize, "%.5f", mview->mol->cell->cellsigma[row]);
3214                                                 break;
3215                                         case 3: case 4: case 5:
3216                                                 if (column == 1)
3217                                                         snprintf(buf, bufsize, "%.4f", mview->mol->cell->cell[row]);
3218                                                 else if (column == 2 && mview->mol->cell->has_sigma)
3219                                                         snprintf(buf, bufsize, "%.4f", mview->mol->cell->cellsigma[row]);
3220                                                 break;
3221                                         case 6: case 7: case 8:
3222                                                 if (column == 1)
3223                                                         snprintf(buf, bufsize, "%d", (int)mview->mol->cell->flags[row - 6]);
3224                                                 break;
3225                                         case 9: case 10: case 11: case 12: {
3226                                                 Vector *vp = (row == 12 ? &mview->mol->cell->origin : &mview->mol->cell->axes[row - 9]);
3227                                                 Double dval = (column == 1 ? vp->x : (column == 2 ? vp->y : vp->z));
3228                                                 snprintf(buf, bufsize, "%.5f", dval);
3229                                                 break;
3230                                         }
3231                                 }
3232                         }
3233                 }
3234         } else if (mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* fractional coords */
3235                 Vector fract_r;
3236                 idx = MainView_tableRowToIndex(mview, row);
3237                 ap[0] = ATOM_AT_INDEX(mol->atoms, idx);
3238                 if (mview->mol->cell != NULL)
3239                         TransformVec(&fract_r, mview->mol->cell->rtr, &(ap[0]->r));
3240                 else fract_r = ap[0]->r;
3241                 switch (column) {
3242                         case 0: snprintf(buf, bufsize, "%d", idx); break;
3243                         case 1: snprintf(buf, bufsize, "%.4s", ap[0]->aname); break;
3244                         case 2: snprintf(buf, bufsize, "%.2s", ap[0]->element); break;
3245                         case 3: snprintf(buf, bufsize, "%.3f", fract_r.x); break;
3246                         case 4: snprintf(buf, bufsize, "%.3f", fract_r.y); break;
3247                         case 5: snprintf(buf, bufsize, "%.3f", fract_r.z); break;
3248                         case 6: snprintf(buf, bufsize, "%.3f", ap[0]->occupancy); break;
3249                         case 7: snprintf(buf, bufsize, "%.3f", ap[0]->tempFactor); break;
3250                         default: buf[0] = 0; break;
3251                 }
3252         } else if (mview->tableIndex == kMainViewMOTableIndex) { /* MO energy */
3253                 BasisSet *bset = mview->mol->bset;
3254                 idx = row;
3255                 buf[0] = 0;
3256                 if (bset != NULL && idx >= 0 && idx < bset->ncomps) {
3257                         if (column == 0) {
3258                                 snprintf(buf, bufsize, "%d", idx + 1);
3259                         } else if (column >= 3) {
3260                                 return;  /*  Cannot happen  */
3261                         } else {
3262                                 /*  column == 1, alpha; column == 2, beta  */
3263                                 char eflag = ' ';
3264                                 if (column == 1) {
3265                                         /*  Alpha  */
3266                                         if (idx >= bset->ne_alpha)
3267                                                 eflag = '*';  /*  Unoccupied  */
3268                                         else {
3269                                                 if (bset->rflag == 2 && idx >= bset->ne_beta)
3270                                                         eflag = 'S';  /*  Singly occupied  */
3271                                         }
3272                                 } else {
3273                                         /*  Beta  */
3274                                         if (bset->rflag != 0)
3275                                                 return;  /*  No beta orbitals  */
3276                                         if (idx >= bset->ne_beta)
3277                                                 eflag = '*';  /*  Unoccupied  */
3278                                         idx += bset->ncomps;
3279                                 }
3280                                 snprintf(buf, bufsize, "%c %.8f", eflag, bset->moenergies[idx]);
3281                         }
3282                 }
3283         }
3284 }
3285
3286 /*  Set color for the table  */
3287 int
3288 MainView_setColorForTable(MainView *mview, int column, int row, float *fg, float *bg)
3289 {
3290         int parType = -1;
3291         int idx;
3292         UnionPar *up;
3293
3294         if (mview->tableIndex == kMainViewParameterTableIndex && column == -1) {
3295                 int src = ParameterTableGetItemSource(mview->mol->par, row);
3296                 if (src == -2) {  /* separator line */
3297                         bg[0] = bg[1] = bg[2] = 0.6;
3298                         return 2;
3299                 } else if (src == -1) { /*  undefined parameters  */
3300                         bg[0] = 1.0;
3301                         bg[1] = bg[2] = 0.2;
3302                         return 2;
3303                 } else if (src == 0) {  /*  local parameter  */
3304                         bg[0] = bg[1] = 1.0;
3305                         bg[2] = 0.6;
3306                         return 2;
3307                 }
3308         } else if (mview->tableIndex == kMainViewMOTableIndex && column == -1) {
3309                 BasisSet *bset = mview->mol->bset;
3310                 int n = 0;
3311                 if (bset == NULL)
3312                         return 0;
3313                 if (row < 0 || row >= bset->ncomps)
3314                         return 0;
3315                 if (row < bset->ne_beta)
3316                         n = 0;  /*  Occupied  */
3317                 else if (row < bset->ne_alpha)
3318                         n = 1;  /*  Half-occupied  */
3319                 else n = 2;  /*  Unoccupied  */
3320                 switch (n) {
3321                         case 1: bg[0] = bg[1] = bg[2] = 0.9; break;
3322                         case 2: bg[0] = bg[1] = bg[2] = 0.8; break;
3323                         default: bg[0] = bg[1] = bg[2] = 1.0; break;
3324                 }
3325                 return 2;
3326         } else if (mview->tableIndex < kMainViewBondTableIndex || mview->tableIndex > kMainViewImproperTableIndex)
3327                 return 0;
3328
3329         /*  Bond etc. table  */
3330         switch (mview->tableIndex) {
3331                 case kMainViewBondTableIndex:
3332                         parType = kBondParType;
3333                         break;
3334                 case kMainViewAngleTableIndex:
3335                         parType = kAngleParType;
3336                         break;
3337                 case kMainViewDihedralTableIndex:
3338                         parType = kDihedralParType;
3339                         break;
3340                 case kMainViewImproperTableIndex:
3341                         parType = kImproperParType;
3342                         break;
3343                 default:
3344                         return 0;
3345         }
3346         if (column < 3)
3347                 return 0;
3348
3349         idx = MainView_tableRowToIndex(mview, row);
3350         up = sParameterOfTypeAtIndex(mview->mol, parType, idx);
3351         if (up == NULL)
3352                 return 0;
3353
3354         if (column == 3) {
3355                 /*  Value column; warn if the value is "abnormal"  */
3356                 int f;
3357                 switch (parType) {
3358                         case kBondParType: f = md_check_abnormal_bond(mview->mol->arena, mview->mol, idx); break;
3359                         case kAngleParType: f = md_check_abnormal_angle(mview->mol->arena, mview->mol, idx); break;
3360                         case kDihedralParType: f = md_check_abnormal_dihedral(mview->mol->arena, mview->mol, idx); break;
3361                         case kImproperParType: f = md_check_abnormal_improper(mview->mol->arena, mview->mol, idx); break;
3362                         default: return 0;
3363                 }
3364                 if (f == 1) {
3365                         bg[0] = 1.0;
3366                         bg[1] = bg[2] = 0.5;
3367                         return 2;
3368                 } else return 0;
3369         } else {
3370                 if (up->bond.src == 0) {
3371                         bg[0] = bg[1] = 1.0;
3372                         bg[2] = 0.6;
3373                         return 2;
3374                 } else if (up->bond.src == -1) {
3375                         bg[0] = 1.0;
3376                         bg[1] = bg[2] = 0.2;
3377                         return 2;
3378                 }
3379                 return 0;
3380         }
3381 }
3382
3383 void
3384 MainView_setSelectionFromTable(MainView *mview)
3385 {
3386         IntGroup *ig, *sel;
3387         if (mview == NULL || mview->mol == NULL)
3388                 return;
3389         switch (mview->tableIndex) {
3390                 case kMainViewAtomTableIndex:
3391                 case kMainViewBondTableIndex:
3392                 case kMainViewAngleTableIndex:
3393                 case kMainViewDihedralTableIndex:
3394                 case kMainViewImproperTableIndex:
3395                 case kMainViewXtalCoordTableIndex:
3396                         break;   /*  Continue  */
3397                 default:
3398                         return;  /*  Do nothing  */
3399         }
3400         ig = MainViewCallback_getTableSelection(mview);
3401         sel = IntGroupNew();
3402         if (ig != NULL) {
3403                 int i, i1, i2;
3404                 Int *ip;
3405                 for (i = 0; (i1 = IntGroupGetNthPoint(ig, i)) >= 0; i++) {
3406                         i2 = IntGroupGetNthPoint(mview->tableCache, i1);
3407                         if (i2 < 0)
3408                                 continue;
3409                         if (mview->tableIndex == kMainViewAtomTableIndex ||
3410                                 mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Atoms */
3411                                 IntGroupAdd(sel, i2, 1);
3412                         } else if (mview->tableIndex == kMainViewBondTableIndex) {  /* Bonds */
3413                                 ip = mview->mol->bonds + i2 * 2;
3414                                 IntGroupAdd(sel, ip[0], 1);
3415                                 IntGroupAdd(sel, ip[1], 1);
3416                         } else if (mview->tableIndex == kMainViewAngleTableIndex) {  /* Angles */
3417                                 ip = mview->mol->angles + i2 * 3;
3418                                 IntGroupAdd(sel, ip[0], 1);
3419                                 IntGroupAdd(sel, ip[1], 1);
3420                                 IntGroupAdd(sel, ip[2], 1);
3421                         } else if (mview->tableIndex == kMainViewDihedralTableIndex || mview->tableIndex == kMainViewImproperTableIndex) {  /* Dihedrals, impropers */
3422                                 ip = (mview->tableIndex == kMainViewDihedralTableIndex ? mview->mol->dihedrals : mview->mol->impropers) + i2 * 4;
3423                                 IntGroupAdd(sel, ip[0], 1);
3424                                 IntGroupAdd(sel, ip[1], 1);
3425                                 IntGroupAdd(sel, ip[2], 1);
3426                                 IntGroupAdd(sel, ip[3], 1);
3427                         }
3428                 }
3429                 IntGroupRelease(ig);
3430         }
3431         MoleculeSetSelection(mview->mol, sel);
3432         IntGroupRelease(sel);
3433 }
3434
3435 static void
3436 sMainView_ParameterTableSetItemText(Molecule *mol, int column, int row, const char *s)
3437 {
3438         Parameter *par;
3439         int type, idx;
3440         const char *kstr = NULL;
3441         UnionPar *up = NULL;
3442         if (mol == NULL || (par = mol->par) == NULL || row < 0)
3443                 return;
3444         up = ParameterTableGetItemPtr(par, row, &type);
3445         if (up == NULL)
3446                 return;
3447         switch (type) {
3448                 case kVdwParType:
3449                         switch (column) {
3450                                 case 1: kstr = "atom_type"; break;
3451                                 case 2: kstr = "eps"; break;
3452                                 case 3: kstr = "r_eq"; break;
3453                                 case 4: kstr = "eps14"; break;
3454                                 case 5: kstr = "r_eq14"; break;
3455                                 case 6: kstr = "atomic_number"; break;
3456                                 case 7: kstr = "weight"; break;
3457                         }
3458                         break;
3459                 case kBondParType:
3460                         switch (column) {
3461                                 case 1: kstr = "atom_types"; break;
3462                                 case 2: kstr = "k"; break;
3463                                 case 3: kstr = "r0"; break;
3464                         }
3465                         break;
3466                 case kAngleParType:
3467                         switch (column) {
3468                                 case 1: kstr = "atom_types"; break;
3469                                 case 2: kstr = "k"; break;
3470                                 case 3: kstr = "a0"; break;
3471                         }
3472                         break;
3473                 case kDihedralParType:
3474                 case kImproperParType:
3475                         switch (column) {
3476                                 case 1: kstr = "atom_types"; break;
3477                                 case 2: kstr = "k"; break;
3478                                 case 3: kstr = "period"; break;
3479                                 case 4: kstr = "phi0"; break;
3480                         }
3481                         break;
3482                 case kVdwPairParType:
3483                         switch (column) {
3484                                 case 1: kstr = "atom_types"; break;
3485                                 case 2: kstr = "eps"; break;
3486                                 case 3: kstr = "r_eq"; break;
3487                                 case 4: kstr = "eps14"; break;
3488                                 case 5: kstr = "r_eq14"; break;
3489                         }
3490                         break;
3491         }
3492         if (column == 9)
3493                 kstr = "comment";
3494         if (kstr == NULL)
3495                 return;
3496         idx = ParameterTableGetItemIndex(par, row, &type);
3497         MolActionCreateAndPerform(mol, SCRIPT_ACTION("iissi"), "set_parameter_attr", type, idx, kstr, s, 0);
3498 }
3499
3500 void
3501 MainView_setValueForTable(MainView *mview, int column, int row, const char *buf)
3502 {
3503         int idx;
3504         char *key;
3505         Molecule *mol;
3506         char temp[256];
3507         if (mview == NULL || (mol = mview->mol) == NULL)
3508                 return;
3509         MainView_valueForTable(mview, column, row, temp, sizeof temp);
3510         if (strcmp(buf, temp) == 0 || buf[0] == 0)
3511                 return;  /*  No change  */
3512         idx = MainView_tableRowToIndex(mview, row);
3513         if (mview->tableIndex == kMainViewAtomTableIndex) { /* Atoms */
3514                 switch (column) {
3515                         case 1: key = "name"; break;
3516                         case 2: key = "atom_type"; break;
3517                         case 3: key = "element"; break;
3518                         case 4: {  /*  Residue  */
3519                                 IntGroup *ig = IntGroupNewWithPoints(idx, 1, -1);
3520                                 MolActionCreateAndPerform(mol, SCRIPT_ACTION("Gs"), "assign_residue", ig, buf);
3521                                 IntGroupRelease(ig);
3522                                 return;
3523                         }
3524                         case 5: key = "x"; break;
3525                         case 6: key = "y"; break;
3526                         case 7: key = "z"; break;
3527                         case 8: key = "charge"; break;
3528                         default: return;
3529                 }
3530                 MolActionCreateAndPerform(mol, SCRIPT_ACTION("iss"), "set_atom_attr", idx, key, buf);
3531         } else if (mview->tableIndex == kMainViewParameterTableIndex) { /* Parameters */
3532                 sMainView_ParameterTableSetItemText(mview->mol, column, row, buf);
3533         } else if (mview->tableIndex == kMainViewUnitCellTableIndex) { /* Unit cell */
3534                 Double cellPars[12], dval;
3535                 char flags[3];
3536                 int n1;
3537                 Vector vecs[4];
3538                 if (mol->cell == NULL) {
3539                         /*  Create a cell  */
3540                         static const Double sCellPars[] = {1, 1, 1, 90, 90, 90};
3541                         MolActionCreateAndPerform(mol, gMolActionSetCell, 6, sCellPars, 0);
3542                 }
3543                 if (row >= 0 && row < 6) {
3544                         n1 = 6;
3545                         memmove(cellPars, mol->cell->cell, sizeof(Double) * 6);
3546                         if (mol->cell->has_sigma)
3547                                 memmove(cellPars + 6, mol->cell->cellsigma, sizeof(Double) * 6);
3548                         else memset(cellPars + 6, 0, sizeof(Double) * 6);
3549                         dval = strtod(buf, NULL);
3550                         if (column == 1) {
3551                                 if (dval == 0.0)
3552                                         return;
3553                                 cellPars[row] = dval;
3554                         } else {
3555                                 n1 = 12;
3556                                 cellPars[row + 6] = dval;
3557                         }
3558                         MolActionCreateAndPerform(mol, gMolActionSetCell, n1, cellPars, 0);
3559                 } else {
3560                         memmove(vecs, mol->cell->axes, sizeof(Vector) * 3);
3561                         vecs[3] = mol->cell->origin;
3562                         memmove(flags, mol->cell->flags, 3);
3563                         if (row >= 6 && row < 9)
3564                                 flags[row - 6] = strtol(buf, NULL, 0);
3565                         else {
3566                                 Vector *vp = vecs + (row - 9);
3567                                 dval = strtod(buf, NULL);
3568                                 if (column == 1)
3569                                         vp->x = dval;
3570                                 else if (column == 2)
3571                                         vp->y = dval;
3572                                 else vp->z = dval;
3573                         }
3574                         MolActionCreateAndPerform(mol, gMolActionSetBox, vecs, vecs + 1, vecs + 2, vecs + 3, (flags[0] != 0) * 4 + (flags[1] != 0) * 2 + (flags[2] != 0), 0);
3575                 }
3576
3577         } else if (mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Fractional coords */
3578                 switch (column) {
3579                         case 1: key = "name"; break;
3580                         case 2: key = "element"; break;
3581                         case 3: key = "fract_x"; break;
3582                         case 4: key = "fract_y"; break;
3583                         case 5: key = "fract_z"; break;
3584                         case 6: key = "occupancy"; break;
3585                         case 7: key = "temp_factor"; break;
3586                         default: return;
3587                 }
3588                 MolActionCreateAndPerform(mol, SCRIPT_ACTION("iss"), "set_atom_attr", idx, key, buf);
3589         }
3590 }
3591
3592 int
3593 MainView_isTableItemEditable(MainView *mview, int column, int row)
3594 {
3595         if (mview == NULL)
3596                 return ParameterTableIsItemEditable(gBuiltinParameters, column, row);
3597         if (mview->mol == NULL)
3598                 return 0;
3599         if (mview->tableIndex == kMainViewParameterTableIndex)
3600                 return ParameterTableIsItemEditable(mview->mol->par, column, row);
3601         if (mview->tableIndex == kMainViewUnitCellTableIndex) {
3602                 if ((row >= 0 && row < 6 && column >= 1 && column <= 3) ||
3603                         (row >= 6 && row < 9 && column == 1) ||
3604                         (row >= 9 && row < 13 && column >= 1 && column <= 3))
3605                         return 1;
3606                 else return 0;
3607         }
3608         if (mview->tableIndex >= 0 && mview->tableIndex < sizeof(sColumnInfo) / sizeof(sColumnInfo[0]))
3609                 return sColumnInfo[mview->tableIndex][column].editable != 0;
3610         else return 0;
3611 }
3612
3613 int
3614 MainView_tableType(MainView *mview)
3615 {
3616         if (mview == NULL)
3617                 return kMainViewParameterTableIndex;
3618         if (mview->mol == NULL)
3619                 return -1;
3620         return mview->tableIndex;
3621 }
3622
3623 void
3624 MainView_dragTableSelectionToRow(MainView *mview, int row)
3625 {
3626         Int *new2old, i, n, count, natoms, start_idx, to_idx;
3627         IntGroup *sel;
3628         if (mview == NULL || mview->mol == NULL)
3629                 return;
3630         if (mview->tableIndex != kMainViewAtomTableIndex && mview->tableIndex != kMainViewXtalCoordTableIndex)
3631                 return;  /*  Only atom tables can respond  */
3632         if (md_is_running(mview->mol->arena)) {
3633                 MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
3634                 return;
3635         }
3636         n = MainView_numberOfRowsInTable(mview);
3637         if (row < 0 || row > n)
3638                 return;  /*  Out of range  */
3639         natoms = mview->mol->natoms;
3640         if (row == n)
3641                 to_idx = natoms;
3642         else
3643                 to_idx = MainView_tableRowToIndex(mview, row);
3644         sel = MoleculeGetSelection(mview->mol);
3645         if (sel == NULL || (count = IntGroupGetCount(sel)) == 0)
3646                 return;
3647         new2old = (Int *)calloc(sizeof(Int), natoms);
3648         if (new2old == NULL)
3649                 return;
3650
3651         //  Place the atoms above the target position
3652         for (i = n = 0; i < to_idx; i++) {
3653                 if (IntGroupLookupPoint(sel, i) < 0)
3654                         new2old[n++] = i;
3655         }
3656         start_idx = n;
3657         //  Place the atoms within the selection
3658         for (i = 0; i < count; i++) {
3659                 new2old[n++] = IntGroupGetNthPoint(sel, i);
3660         }
3661         //  Place the remaining atoms
3662         for (i = to_idx; i < natoms; i++) {
3663                 if (IntGroupLookupPoint(sel, i) < 0)
3664                         new2old[n++] = i;
3665         }
3666         MolActionCreateAndPerform(mview->mol, gMolActionRenumberAtoms, n, new2old);
3667
3668         //  Change selection
3669         sel = IntGroupNewWithPoints(start_idx, count, -1);
3670         MolActionCreateAndPerform(mview->mol, gMolActionSetSelection, sel);
3671         IntGroupRelease(sel);
3672 }
3673
3674 IntGroup *
3675 MainView_selectedMO(MainView *mview)
3676 {
3677         if (mview == NULL || mview->mol == NULL || mview->tableIndex != kMainViewMOTableIndex)
3678                 return NULL;
3679         return MainViewCallback_getTableSelection(mview);  /*  Note: the indices are 0 based  */
3680 }