MolLib/MainView.c

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