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