MolLib/MainView.c

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