OSDN Git Service

MO surface display is improved
[molby/Molby.git] / MolLib / MainView.c
index 4f4d953..78d360f 100755 (executable)
  GNU General Public License for more details.
  */
 
+/*  On MinGW, #include of <GL/gl.h> should be before something in MolLib.h.
+ Therefore, we include MainView.h first separately and then include
+ the rest of MolLib.h.  */
+#include "MainView.h"
 #include "MolLib.h"
+#include "Ruby_bind/Molby_extern.h"
 
 #include "MD/MDCore.h"
 #include "MD/MDGraphite.h"
 #include <stdlib.h>
 #include <string.h>
 #include <ctype.h>
+#include <math.h>
 
 #define biso2radius(r) ((r) > 0.5 ? sqrt((r) / 78.9568352087147) : 0.08)
 
+/*  pp is a pointer to Parameter record  */
+#define VDW_RADIUS(pp) (pp->vdw_radius == 0.0 ? pp->radius * 0.51 + 1.20 : pp->vdw_radius)
+
 /*  Invalid bond/angle/torsion value, used in internal cache  */
 const Double kInvalidValue = -10000000.0;
 
-#pragma mark ==== MainView public methods ====
+#pragma mark ==== OpenGL utility functions ===
 
-MainView *
-MainView_newMainView(void *ref)
+static void __gluMultMatrixVecf(const GLdouble matrix[16], const GLdouble in[4], GLdouble out[4])
 {
-       static GLdouble sIdentity[16] = {
-               1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1
-       };
-       MainView *mview = (MainView *)malloc(sizeof(MainView));
-       if (mview != NULL) {
-               memset(mview, 0, sizeof(MainView));
-               mview->ref = ref;
-               mview->track = TrackballNew();
-               mview->mode = kTrackballRotateMode;
-               mview->tempAtoms[0] = mview->tempAtoms[1] = -1;
-               memmove(mview->modelview_matrix, sIdentity, sizeof(sIdentity));
-               memmove(mview->projection_matrix, sIdentity, sizeof(sIdentity));
-               mview->atomRadius = 0.4;
-               mview->bondRadius = 0.1;
-               mview->probabilityScale = 1.5382;
-               mview->dimension = 10.0;
-       /*      MainView_resizeToFit(mview); */
-               mview->showHydrogens = mview->showDummyAtoms = mview->showExpandedAtoms = 1;
-               mview->showPeriodicBox = 1;
-               mview->showGraphite = 5;
-               mview->tableCache = IntGroupNew();
-               mview->tableSelection = IntGroupNew();
-       }
-       return mview;
+    int i;
+    for (i = 0; i < 4; i++) {
+        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];
+    }
 }
 
-void
-MainView_release(MainView *mview)
+static int __gluInvertMatrixf(const GLdouble m[16], GLdouble invOut[16])
 {
-       int i;
-       if (mview != NULL) {
-               MainView_setMolecule(mview, NULL);
-               TrackballRelease(mview->track);
-               IntGroupRelease(mview->tableCache);
-               IntGroupRelease(mview->tableSelection);
-               if (mview->nlabels > 0) {
-                       for (i = 0; i < mview->nlabels; i++) {
-                               MainViewCallback_releaseLabel(mview->labels[i].label);
-                       }
-                       free(mview->labels);
-                       free(mview->sortedLabels);
-               }
-               if (mview->rotateFragment != NULL)
-                       IntGroupRelease(mview->rotateFragment);
-               if (mview->rotateFragmentOldPos != NULL)
-                       free(mview->rotateFragmentOldPos);
-               if (mview->visibleFlags != NULL)
-                       free(mview->visibleFlags);
-       }
-       free(mview);
+    GLdouble inv[16], det;
+    int i;
+    
+    inv[0] =   m[5]*m[10]*m[15] - m[5]*m[11]*m[14] - m[9]*m[6]*m[15]
+    + m[9]*m[7]*m[14] + m[13]*m[6]*m[11] - m[13]*m[7]*m[10];
+    inv[4] =  -m[4]*m[10]*m[15] + m[4]*m[11]*m[14] + m[8]*m[6]*m[15]
+    - m[8]*m[7]*m[14] - m[12]*m[6]*m[11] + m[12]*m[7]*m[10];
+    inv[8] =   m[4]*m[9]*m[15] - m[4]*m[11]*m[13] - m[8]*m[5]*m[15]
+    + m[8]*m[7]*m[13] + m[12]*m[5]*m[11] - m[12]*m[7]*m[9];
+    inv[12] = -m[4]*m[9]*m[14] + m[4]*m[10]*m[13] + m[8]*m[5]*m[14]
+    - m[8]*m[6]*m[13] - m[12]*m[5]*m[10] + m[12]*m[6]*m[9];
+    inv[1] =  -m[1]*m[10]*m[15] + m[1]*m[11]*m[14] + m[9]*m[2]*m[15]
+    - m[9]*m[3]*m[14] - m[13]*m[2]*m[11] + m[13]*m[3]*m[10];
+    inv[5] =   m[0]*m[10]*m[15] - m[0]*m[11]*m[14] - m[8]*m[2]*m[15]
+    + m[8]*m[3]*m[14] + m[12]*m[2]*m[11] - m[12]*m[3]*m[10];
+    inv[9] =  -m[0]*m[9]*m[15] + m[0]*m[11]*m[13] + m[8]*m[1]*m[15]
+    - m[8]*m[3]*m[13] - m[12]*m[1]*m[11] + m[12]*m[3]*m[9];
+    inv[13] =  m[0]*m[9]*m[14] - m[0]*m[10]*m[13] - m[8]*m[1]*m[14]
+    + m[8]*m[2]*m[13] + m[12]*m[1]*m[10] - m[12]*m[2]*m[9];
+    inv[2] =   m[1]*m[6]*m[15] - m[1]*m[7]*m[14] - m[5]*m[2]*m[15]
+    + m[5]*m[3]*m[14] + m[13]*m[2]*m[7] - m[13]*m[3]*m[6];
+    inv[6] =  -m[0]*m[6]*m[15] + m[0]*m[7]*m[14] + m[4]*m[2]*m[15]
+    - m[4]*m[3]*m[14] - m[12]*m[2]*m[7] + m[12]*m[3]*m[6];
+    inv[10] =  m[0]*m[5]*m[15] - m[0]*m[7]*m[13] - m[4]*m[1]*m[15]
+    + m[4]*m[3]*m[13] + m[12]*m[1]*m[7] - m[12]*m[3]*m[5];
+    inv[14] = -m[0]*m[5]*m[14] + m[0]*m[6]*m[13] + m[4]*m[1]*m[14]
+    - m[4]*m[2]*m[13] - m[12]*m[1]*m[6] + m[12]*m[2]*m[5];
+    inv[3] =  -m[1]*m[6]*m[11] + m[1]*m[7]*m[10] + m[5]*m[2]*m[11]
+    - m[5]*m[3]*m[10] - m[9]*m[2]*m[7] + m[9]*m[3]*m[6];
+    inv[7] =   m[0]*m[6]*m[11] - m[0]*m[7]*m[10] - m[4]*m[2]*m[11]
+    + m[4]*m[3]*m[10] + m[8]*m[2]*m[7] - m[8]*m[3]*m[6];
+    inv[11] = -m[0]*m[5]*m[11] + m[0]*m[7]*m[9] + m[4]*m[1]*m[11]
+    - m[4]*m[3]*m[9] - m[8]*m[1]*m[7] + m[8]*m[3]*m[5];
+    inv[15] =  m[0]*m[5]*m[10] - m[0]*m[6]*m[9] - m[4]*m[1]*m[10]
+    + m[4]*m[2]*m[9] + m[8]*m[1]*m[6] - m[8]*m[2]*m[5];
+    
+    det = m[0]*inv[0] + m[1]*inv[4] + m[2]*inv[8] + m[3]*inv[12];
+    if (det == 0)
+        return GL_FALSE;
+    
+    det = 1.0 / det;
+    
+    for (i = 0; i < 16; i++)
+        invOut[i] = inv[i] * det;
+    
+    return GL_TRUE;
 }
 
+static void __gluMultMatricesf(const GLdouble a[16], const GLdouble b[16],
+                               GLdouble r[16])
+{
+    int i, j;
+    for (i = 0; i < 4; i++) {
+        for (j = 0; j < 4; j++) {
+            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];
+        }
+    }
+}
+
+GLint
+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)
+{
+    GLdouble in[4];
+    GLdouble out[4];
+    
+    in[0] = objx;
+    in[1] = objy;
+    in[2] = objz;
+    in[3] = 1.0;
+    __gluMultMatrixVecf(modelMatrix, in, out);
+    __gluMultMatrixVecf(projMatrix, out, in);
+    if (in[3] == 0.0)
+        return(GL_FALSE);
+    in[0] /= in[3];
+    in[1] /= in[3];
+    in[2] /= in[3];
+    /* Map x, y and z to range 0-1 */
+    in[0] = in[0] * 0.5 + 0.5;
+    in[1] = in[1] * 0.5 + 0.5;
+    in[2] = in[2] * 0.5 + 0.5;
+    /* Map x,y to viewport */
+    in[0] = in[0] * viewport[2] + viewport[0];
+    in[1] = in[1] * viewport[3] + viewport[1];
+    *winx = in[0];
+    *winy = in[1];
+    *winz = in[2];
+    return(GL_TRUE);
+}
+
+GLint
+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)
+{
+    GLdouble finalMatrix[16];
+    GLdouble in[4];
+    GLdouble out[4];
+    
+    __gluMultMatricesf(modelMatrix, projMatrix, finalMatrix);
+    if (!__gluInvertMatrixf(finalMatrix, finalMatrix))
+        return(GL_FALSE);
+    
+    in[0] = winx;
+    in[1] = winy;
+    in[2] = winz;
+    in[3] = 1.0;
+    
+    /* Map x and y from window coordinates */
+    in[0] = (in[0] - viewport[0]) / viewport[2];
+    in[1] = (in[1] - viewport[1]) / viewport[3];
+    
+    /* Map to range -1 to 1 */
+    in[0] = in[0] * 2 - 1;
+    in[1] = in[1] * 2 - 1;
+    in[2] = in[2] * 2 - 1;
+    
+    __gluMultMatrixVecf(finalMatrix, in, out);
+    if (out[3] == 0.0) return(GL_FALSE);
+    out[0] /= out[3];
+    out[1] /= out[3];
+    out[2] /= out[3];
+    *objx = out[0];
+    *objy = out[1];
+    *objz = out[2];
+    return(GL_TRUE);
+}
+
+void
+myGluPerspective(GLdouble fovy_degree, GLdouble aspect, GLdouble zNear, GLdouble zFar) {
+    GLdouble m[16];
+    GLdouble fovy_rad = fovy_degree * (3.14159265358979 / 180.0);
+    GLdouble f = 1.0 / tan(fovy_rad / (GLdouble)2.0);
+    m[0] = f / aspect;
+    m[1] = 0.0;
+    m[2] = 0.0;
+    m[3] = 0.0;
+    m[4] = 0.0;
+    m[5] = f;
+    m[6] = 0.0;
+    m[7] = 0.0;
+    m[8] = 0.0;
+    m[9] = 0.0;
+    m[10] = (zFar + zNear) / (zNear - zFar);
+    m[11] = -1.0;
+    m[12] = 0.0;
+    m[13] = 0.0;
+    m[14] = (2 * zFar * zNear) / (zNear - zFar);
+    m[15] = 0.0;
+    glMatrixMode(GL_PROJECTION);
+    glLoadMatrixd(m);
+}
+
+#pragma mark ==== MainView public methods ====
+
 void
-MainView_setMolecule(MainView *mview, struct Molecule *mol)
+MainView_setViewObject(MainView *mview, void *ref)
 {
-       if (mview == NULL || mview->mol == mol)
-               return;
-       if (mview->mol != NULL) {
-               mview->mol->mview = NULL;  /*  No need to release  */
-               MoleculeRelease(mview->mol);
-       }
-       mview->mol = mol;
-       if (mol != NULL) {
-               MoleculeRetain(mol);
-               mol->mview = mview;  /*  No retain  */
-               MainViewCallback_moleculeReplaced(mview, mol);
-               MainView_resizeToFit(mview);
-               MoleculeCallback_notifyModification(mol, 0);
-       /*      MainViewCallback_setNeedsDisplay(mview, 1); */
+       static GLdouble sIdentity[16] = {
+               1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1
+       };
+       if (mview != NULL) {
+               if (mview->ref == NULL) {
+                       /*  Initialize GUI-only members  */
+                       mview->mode = kTrackballRotateMode;
+                       mview->tempAtoms[0] = mview->tempAtoms[1] = -1;
+                       memmove(mview->modelview_matrix, sIdentity, sizeof(sIdentity));
+                       memmove(mview->projection_matrix, sIdentity, sizeof(sIdentity));
+                       mview->tableCache = IntGroupNew();
+                       mview->tableSelection = IntGroupNew();
+               }
+               mview->ref = ref;  /*  No retain  */
+               IntGroupClear(mview->tableCache);
+               IntGroupClear(mview->tableSelection);
+        if (mview->ref != NULL)
+            MoleculeCallback_notifyModification(mview->mol, 0);
        }
 }
 
@@ -177,7 +297,8 @@ MainView_refreshCachedInfo(MainView *mview)
        IntGroupClear(mview->tableCache);
        IntGroupClear(mview->tableSelection);
                
-       if (mview->tableIndex == kMainViewAtomTableIndex) {  /* Atoms */
+       if (mview->tableIndex == kMainViewAtomTableIndex ||
+               mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Atoms */
                for (i = j = 0; i < mol->natoms; i++) {
                        f1 = mview->visibleFlags[i];
                        if ((f1 & 1) != 0) {
@@ -249,18 +370,26 @@ MainView_refreshCachedInfo(MainView *mview)
                                j++;
                        }
                }
-       } else if (mview->tableIndex == kMainViewParameterTableIndex) {  /* Parameter infos */
-               /*  Do nothing (tableCache will not be used)  */
-       } else if (mview->tableIndex == kMainViewMOTableIndex) {  /* MO infos  */
-               /*  Really no need to cache info, but create it anyway to simplify code  */
-               if (mol->bset != NULL && mol->bset->nmos > 0)
-                       IntGroupAdd(mview->tableCache, 0, mol->bset->nmos);
+       } else {
+               /*  Cache is left empty (not used)  */
        }
+
+//     } else if (mview->tableIndex == kMainViewParameterTableIndex ||  /* Parameter infos */
+//                        mview->tableIndex == kMainViewUnitCellTableIndex) {   /* Unit cell infos */
+//             /*  Do nothing (tableCache will not be used)  */
+//     } else if (mview->tableIndex == kMainViewMOTableIndex) {  /* MO infos  */
+//             /*  Really no need to cache info, but create it anyway to simplify code  */
+//             if (mol->bset != NULL && mol->bset->ncomps > 0)
+//                     IntGroupAdd(mview->tableCache, 0, mol->bset->ncomps);
+//     }
        
        /*  Clear internal selection flag  */
        for (i = 0; i < mol->natoms; i++) {
                mview->visibleFlags[i] &= ~2;
        }
+    
+    /*  Store the tableIndex value  */
+    mview->cachedTableIndex = mview->tableIndex;
 }
 
 #pragma mark ====== 2D/3D transform operations ======
@@ -269,7 +398,9 @@ void
 MainView_resizeToFit(MainView *mview)
 {
        Vector p;
-       float f[4];
+       double f[4];
+    if (!gUseGUI)
+        return;
        if (mview == NULL || mview->mol == NULL)
                return;
        if (mview->mol->natoms == 0) {
@@ -277,8 +408,6 @@ MainView_resizeToFit(MainView *mview)
                return;
        }
        MoleculeCenterOfMass(mview->mol, &p, NULL);
-/*     if (mview->mol->is_xtal_coord)
-               TransformVec(&p, mview->mol->cell->tr, &p); */
        f[0] = -p.x / mview->dimension;
        f[1] = -p.y / mview->dimension;
        f[2] = -p.z / mview->dimension;
@@ -296,8 +425,6 @@ MainView_resizeToFit(MainView *mview)
                double r0 = 0.0, r1, scale;
                for (i = 0, ap = mview->mol->atoms; i < mview->mol->natoms; i++, ap = ATOM_NEXT(ap)) {
                        q = ap->r;
-               /*      if (mview->mol->is_xtal_coord)
-                               TransformVec(&q, mview->mol->cell->tr, &q); */
                        VecDec(q, p);
                        r1 = VecLength(q);
                        if (r1 > r0)
@@ -314,113 +441,28 @@ MainView_resizeToFit(MainView *mview)
        }
 
        MainViewCallback_setNeedsDisplay(mview, 1);
-
-#if 0
-       GLdouble mm[4][4];
-       GLdouble p0[4];
-       int i, j, natoms;
-       const Atom *ap;
-       const ElementPar *app;
-       GLdouble *pp, *p1;
-       GLdouble max[3], min[3], center[3];
-       float frame[4], width, height, cot15, rsin15, cot_th, rsin_th, d, trans[3];
-       Molecule *mol;
-       
-       if (mview == NULL || mview->mol == NULL)
-               return;
-       mol = mview->mol;
-       
-       /*  Transform the coordinates of all atoms with the rotation part of
-        *  the model-view matrix  */
-       memmove((GLdouble *)mm, mview->modelview_matrix, sizeof(float) * 16);
-       natoms = mol->natoms;
-       pp = (GLdouble *)calloc(sizeof(GLdouble), 4 * natoms);
-       MALLOC_CHECK(pp, "resizing the model to fit");
-       for (i = 0; i < natoms; i++) {
-               ap = ATOM_AT_INDEX(mol->atoms, i);
-               if (ap == NULL)
-                       continue;
-               p0[0] = ap->r.x; p0[1] = ap->r.y; p0[2] = ap->r.z;
-               p0[3] = 0.0;  /*  This disables the translation part  */
-               p1 = &pp[i * 4];
-               MAT_DOT_VEC_4X4(p1, mm, p0);  /* #define'd in GLUT/vvector.h */
-       }
-       
-       /*  Determine the center for each axis  */
-       max[0] = max[1] = max[2] = -1e20;
-       min[0] = min[1] = min[2] = 1e20;
-       for (i = 0; i < natoms; i++) {
-               p1 = &pp[i * 4];
-               for (j = 0; j < 3; j++) {
-                       if (p1[j] > max[j])
-                               max[j] = p1[j];
-                       if (p1[j] < min[j])
-                               min[j] = p1[j];
-               }
-       }
-       for (j = 0; j < 3; j++)
-               center[j] = (min[j] + max[j]) * 0.5;
-       
-       /*  Get the frame size  */
-       MainViewCallback_frame(mview, frame);
-       width = frame[2] - frame[0];
-       height = frame[3] - frame[1];
-
-       /*  Calculate the minimum distance from which the view angle of
-        *  the atom becomes no more than 15 degree (for the y direction)
-        *  and theta degree (for the x direction, theta = arctan(h/w*tan(15 deg))) */
-       cot15 = 3.73205080756888;
-       rsin15 = 3.86370330515628;
-       cot_th = width / height * cot15;
-       rsin_th = sqrt(1.0 + cot_th * cot_th);
-       d = 0;
-       for (i = 0; i < natoms; i++) {
-               float dx, dy, z, r;
-               ap = ATOM_AT_INDEX(mol->atoms, i);
-               if (ap == NULL)
-                       continue;
-               app = &(gBuiltinParameters->atomPars[ap->atomicNumber]);
-               if (app == NULL)
-                       continue;
-               r = app->radius * mview->atomRadius;
-               z = pp[i * 4 + 2] - center[2];
-               dx = pp[i * 4] * cot_th + r * rsin_th + z;
-               dy = pp[i * 4 + 1] * cot15 + r * rsin15 + z;
-               if (d < dx)
-                       d = dx;
-               if (d < dy)
-                       d = dy;
-       }
-       mview->dimension = d / cot15;
-       TrackballSetScale(mview->track, 1.0);
-       trans[0] = center[0] / mview->dimension;
-       trans[1] = center[1] / mview->dimension;
-       trans[2] = center[2] / mview->dimension;
-       TrackballSetTranslate(mview->track, trans);
-
-       free(pp);
-#endif
 }
 
 int
-MainView_convertScreenPositionToObjectPosition(MainView *mview, const GLfloat *screenPos, GLfloat *objectPos)
+MainView_convertScreenPositionToObjectPosition(MainView *mview, const GLfloat *screenPos, double *objectPos)
 {
        float rect[4];
     GLint viewport[4], n;
     GLfloat winZ;
     GLdouble posX, posY, posZ;
+    float scale = mview->view_scale;
        if (mview == NULL)
                return 0;
        MainViewCallback_frame(mview, rect);
     viewport[0] = viewport[1] = 0;
-    viewport[2] = (GLint)(rect[2] - rect[0]);
-    viewport[3] = (GLint)(rect[3] - rect[1]);
+    viewport[2] = (GLint)(rect[2] - rect[0]) * scale;
+    viewport[3] = (GLint)(rect[3] - rect[1]) * scale;
        MainViewCallback_lockFocus(mview);
     if (screenPos[2] >= 0.0)
         winZ = screenPos[2];
     else
-        glReadPixels(screenPos[0], screenPos[1], 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
-    gluUnProject(screenPos[0], screenPos[1], winZ, mview->modelview_matrix, mview->projection_matrix, viewport, &posX, &posY, &posZ);
+        glReadPixels(screenPos[0] * scale, screenPos[1] * scale, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &winZ);
+    myGluUnProject(screenPos[0] * scale, screenPos[1] * scale, winZ, mview->modelview_matrix, mview->projection_matrix, viewport, &posX, &posY, &posZ);
     n = glGetError();
        MainViewCallback_unlockFocus(mview);
        objectPos[0] = posX;
@@ -433,21 +475,22 @@ MainView_convertScreenPositionToObjectPosition(MainView *mview, const GLfloat *s
 }
 
 int
-MainView_convertObjectPositionToScreenPosition(MainView *mview, const GLfloat *objectPos, GLfloat *screenPos)
+MainView_convertObjectPositionToScreenPosition(MainView *mview, const double *objectPos, GLfloat *screenPos)
 {
        float rect[4];
     GLint viewport[4];
     GLdouble objX, objY, objZ;
+    float scale = mview->view_scale;
        if (mview == NULL)
                return 0;
        MainViewCallback_frame(mview, rect);
     viewport[0] = viewport[1] = 0;
-    viewport[2] = (GLint)(rect[2] - rect[0]);
-    viewport[3] = (GLint)(rect[3] - rect[1]);
-    gluProject(objectPos[0], objectPos[1], objectPos[2], mview->modelview_matrix, mview->projection_matrix, viewport, &objX, &objY, &objZ);
+    viewport[2] = (GLint)(rect[2] - rect[0]) * scale;
+    viewport[3] = (GLint)(rect[3] - rect[1]) * scale;
+    myGluProject(objectPos[0], objectPos[1], objectPos[2], mview->modelview_matrix, mview->projection_matrix, viewport, &objX, &objY, &objZ);
     if (glGetError() == GL_NO_ERROR) {
-               screenPos[0] = objX;
-               screenPos[1] = objY;
+               screenPos[0] = objX / scale;
+               screenPos[1] = objY / scale;
                screenPos[2] = objZ;
        /*      fprintf(stderr, "object(%.3f,%.3f,%.3f) screen(%.3f,%.3f,%.3f)\n", objectPos[0], objectPos[1], objectPos[2], screenPos[0], screenPos[1], screenPos[2]); */
                return 1;
@@ -457,7 +500,8 @@ MainView_convertObjectPositionToScreenPosition(MainView *mview, const GLfloat *o
 int
 MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex1, int *outIndex2, int mouseCheck, int ignoreExpandedAtoms)
 {
-       float screenPos[3], op[3], oq[3], pqlen, pqlen2;
+       float screenPos[3];
+       double op[3], oq[3], pqlen, pqlen2;
        Vector pq, pa, v1, r1, r2;
     int i, natoms, nbonds;
        float r, d2, z;
@@ -474,7 +518,6 @@ MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex
        }
        mol = mview->mol;
 
-#if 1
        screenPos[0] = mousePos[0];
        screenPos[1] = mousePos[1];
        screenPos[2] = -1.0;
@@ -496,7 +539,7 @@ MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex
        minDepth = 100.0;
        for (i = 0; i < natoms; i++) {
                Vector pq1, pa1;
-               float pq1len2, pq1len;
+               double pq1len2, pq1len;
                if (mouseCheck && i % 50 == 0 && MainViewCallback_mouseCheck(mview))
                        return 0;  /*  If mouse event is detected return immediately  */
                /*  Examine if an atom is visible or not  */
@@ -516,7 +559,7 @@ MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex
                pa.x = r1.x - op[0];
                pa.y = r1.y - op[1];
                pa.z = r1.z - op[2];
-               if (ap->aniso != NULL) {
+               if (mview->showEllipsoids && ap->aniso != NULL) {
                        /*  Convert to ellipsoid principal axes  */
                        Mat33 m1;
                        Aniso *anp = ap->aniso;
@@ -533,7 +576,7 @@ MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex
                                ep = &(gElementParameters[ap->atomicNumber]);
                                if (ep == NULL)
                                        continue;
-                               r = ep->radius * mview->atomRadius;
+                               r = VDW_RADIUS(ep) * mview->atomRadius;
                        }
                        pa1 = pa;
                        pq1 = pq;
@@ -620,76 +663,6 @@ MainView_findObjectAtPoint(MainView *mview, const float *mousePos, int *outIndex
        *outIndex1 = n1;
        *outIndex2 = n2;
        return (n1 >= 0 || n2 >= 0);
-               
-#else
-       /*  Examine if anything is drawn at the point  */
-       screenPos[0] = mousePos[0];
-       screenPos[1] = mousePos[1];
-       screenPos[2] = -1.0;
-       if (MainView_convertScreenPositionToObjectPosition(mview, screenPos, objectPos) == 0)
-               return 0;  /*  Nothing is here  */
-       op[0] = objectPos[0];
-       op[1] = objectPos[1];
-       op[2] = objectPos[2];
-
-       /*  Examine the distance from the atom center with op, and select the one
-           with the smallest difference with the radius on the screen  */
-       atomRadius = mview->atomRadius;
-       bondRadius = mview->bondRadius;
-       natoms = mol->natoms;
-       mindr = 100.0;
-       n1 = n2 = -1;
-    for (i = 0; i < natoms; i++) {
-        ap = ATOM_AT_INDEX(mol->atoms, i);
-               if (ap == NULL)
-                       continue;
-               dp = &(gBuiltinParameters->atomPars[ap->atomicNumber]);
-               if (dp == NULL)
-                       continue;
-               p[0] = objectPos[0] - ap->r.x;
-               p[1] = objectPos[1] - ap->r.y;
-               p[2] = objectPos[2] - ap->r.z;
-               rr = sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]) / (dp->radius * atomRadius);
-               dr = fabs(1.0 - rr);
-               if (dr < mindr) {
-                       mindr = dr;
-                       n1 = i;
-               }
-       }
-       if (mindr > 0.1)
-               n1 = -1;
-       nbonds = mol->nbonds;
-    for (i = 0; i < nbonds; i++) {
-               float q;
-               ip = &(mol->bonds[i * 2]);
-               ap = ATOM_AT_INDEX(mol->atoms, ip[0]);
-               bp = ATOM_AT_INDEX(mol->atoms, ip[1]);
-               p[0] = objectPos[0] - ap->r.x;
-               p[1] = objectPos[1] - ap->r.y;
-               p[2] = objectPos[2] - ap->r.z;
-               b[0] = bp->r.x - ap->r.x;
-               b[1] = bp->r.y - ap->r.y;
-               b[2] = bp->r.z - ap->r.z;
-               /*  ap.r + q*(b[0],b[1],b[2]) is the closest point on line AB to objectPos */
-               /*  If q is outside [0, 1] then the point is outside the line segment AB
-                   so it is ignored  */
-               q = (p[0]*b[0]+p[1]*b[1]+p[2]*b[2])/(b[0]*b[0]+b[1]*b[1]+b[2]*b[2]);
-               if (q < 0 || q > 1)
-                       continue;
-               p[0] -= q * b[0];
-               p[1] -= q * b[1];
-               p[2] -= q * b[2];
-               dr = fabs(bondRadius - sqrt(p[0]*p[0] + p[1]*p[1] + p[2]*p[2]));
-               if (dr < mindr) {
-                       mindr = dr;
-                       n1 = ip[0];
-                       n2 = ip[1];
-               }
-    }
-       *outIndex1 = n1;
-       *outIndex2 = n2;
-       return (n1 >= 0 || n2 >= 0);
-#endif
 }
 
 int
@@ -699,7 +672,7 @@ MainView_screenCenterPointOfAtom(MainView *mview, int index, float *outScreenPos
        const ElementPar *dp;
        Vector cv, pv, v;
        Double rad, w;
-       float p[3];
+       double p[3];
 
        if (mview == NULL || mview->mol == NULL || index < 0 || index >= mview->mol->natoms)
                return 0;
@@ -742,7 +715,7 @@ MainView_screenCenterPointOfAtom(MainView *mview, int index, float *outScreenPos
                }
        } else {
                dp = &(gElementParameters[ap->atomicNumber]);
-               rad = dp->radius * mview->atomRadius;
+               rad = VDW_RADIUS(dp) * mview->atomRadius;
                w = rad / VecLength(pv) * 1.1;
                VecScaleSelf(pv, w);
        }
@@ -758,6 +731,8 @@ MainView_screenCenterPointOfAtom(MainView *mview, int index, float *outScreenPos
 void
 MainView_getCamera(MainView *mview, Vector *outCamera, Vector *outLookAt, Vector *outUp)
 {
+    if (!gUseGUI)
+        return;
        if (mview != NULL) {
                *outCamera = mview->camera;
                *outLookAt = mview->lookat;
@@ -767,8 +742,24 @@ MainView_getCamera(MainView *mview, Vector *outCamera, Vector *outLookAt, Vector
 
 #pragma mark ====== Draw model ======
 
+static void
+enableLighting(void)
+{
+       static GLfloat pos[] = {0.0f, 0.0f, 1.0f, 0.0f};
+       glEnable(GL_LIGHTING);
+#if __WXMAC__
+       /*  On Mac OS 10.6, redefining the light position seems to be necessary
+               every time lighting is made enabled  */
+       glMatrixMode(GL_MODELVIEW);
+       glPushMatrix();
+       glLoadIdentity();
+       glLightfv(GL_LIGHT0, GL_POSITION, pos);
+       glPopMatrix();
+#endif
+}
+
 void
-MainView_initializeOpenGLView(MainView *mview)
+MainView_initializeOpenGL(void)
 {
        static GLfloat ambient[] = {0.6, 0.6, 0.6, 1.0};  // Some white ambient light.
        static GLfloat diffuse[] = {1.0, 1.0, 1.0, 1.0};  // A white light.
@@ -786,9 +777,6 @@ MainView_initializeOpenGLView(MainView *mview)
        //  Enable blending
        glEnable(GL_BLEND);
        glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
-
-       if (mview != NULL)
-               mview->isInitialized = 1;
 }
 
 /*  Get orthogonal unit vectors  */
@@ -853,10 +841,11 @@ setSinCache(int sect)
 }
 
 static void
-drawCylinder(const GLfloat *a, const GLfloat *b, GLfloat r, int sect)
+drawCylinder(const GLfloat *a, const GLfloat *b, GLfloat r, int sect, int closed)
 {
     GLfloat *c, *s;
     int n, i;
+       float nx, ny, nz;
     GLfloat d[3], v[3], w[3];
     n = setSinCache(sect);
     if (n <= 0)
@@ -870,7 +859,6 @@ drawCylinder(const GLfloat *a, const GLfloat *b, GLfloat r, int sect)
         return;
     glBegin(GL_QUAD_STRIP);
     for (i = 0; i <= n; i++) {
-        float nx, ny, nz;
         nx = v[0] * c[i] + w[0] * s[i];
         ny = v[1] * c[i] + w[1] * s[i];
         nz = v[2] * c[i] + w[2] * s[i];
@@ -878,7 +866,76 @@ drawCylinder(const GLfloat *a, const GLfloat *b, GLfloat r, int sect)
         glVertex3f(a[0] + r * nx, a[1] + r * ny, a[2] + r * nz);
         glVertex3f(b[0] + r * nx, b[1] + r * ny, b[2] + r * nz);
     }
-    glEnd();    
+    glEnd();
+       if (closed) {
+               glBegin(GL_TRIANGLE_FAN);
+               glNormal3f(-d[0], -d[1], -d[2]);
+               for (i = 0; i <= n; i++) {
+                       nx = v[0] * c[i] + w[0] * s[i];
+                       ny = v[1] * c[i] + w[1] * s[i];
+                       nz = v[2] * c[i] + w[2] * s[i];
+                       glVertex3f(a[0] + r * nx, a[1] + r * ny, a[2] + r * nz);
+               }
+               glEnd();
+               glBegin(GL_TRIANGLE_FAN);
+               glNormal3f(d[0], d[1], d[2]);
+               for (i = 0; i <= n; i++) {
+                       nx = v[0] * c[i] + w[0] * s[i];
+                       ny = v[1] * c[i] + w[1] * s[i];
+                       nz = v[2] * c[i] + w[2] * s[i];
+                       glVertex3f(b[0] + r * nx, b[1] + r * ny, b[2] + r * nz);
+               }
+               glEnd();
+       }
+}
+
+static void
+drawCone(const GLfloat *a, const GLfloat *b, GLfloat r, int sect, int closed)
+{
+    GLfloat *c, *s;
+    int n, i;
+    GLfloat d[3], v[3], w[3];
+       Vector p1, nv;
+    n = setSinCache(sect);
+    if (n <= 0)
+        return;
+    s = sSinCache;
+    c = &sSinCache[n/4];
+    d[0] = b[0] - a[0];
+    d[1] = b[1] - a[1];
+    d[2] = b[2] - a[2];
+    if (getOrthogonalVectors(d, v, w) == 0)
+        return;
+    glBegin(GL_TRIANGLE_FAN);
+       nv.x = d[0];
+       nv.y = d[1];
+       nv.z = d[2];
+       NormalizeVec(&nv, &nv);
+       glNormal3f(nv.x, nv.y, nv.z);
+       glVertex3f(b[0], b[1], b[2]);
+    for (i = 0; i <= n; i++) {
+        nv.x = v[0] * c[i] + w[0] * s[i];
+        nv.y = v[1] * c[i] + w[1] * s[i];
+        nv.z = v[2] * c[i] + w[2] * s[i];
+               glNormal3f(nv.x, nv.y, nv.z);
+               p1.x = a[0] + r * nv.x;
+               p1.y = a[1] + r * nv.y;
+               p1.z = a[2] + r * nv.z;
+        glNormal3f(nv.x, nv.y, nv.z);
+               glVertex3f(p1.x, p1.y, p1.z);
+    }
+    glEnd();
+       if (closed) {
+               glBegin(GL_TRIANGLE_FAN);
+               glNormal3f(d[0], d[1], d[2]);
+               for (i = 0; i <= n; i++) {
+                       p1.x = a[0] + r * (v[0] * c[i] + w[0] * s[i]);
+                       p1.y = a[1] + r * (v[1] * c[i] + w[1] * s[i]);
+                       p1.z = a[2] + r * (v[2] * c[i] + w[2] * s[i]);
+                       glVertex3f(p1.x, p1.y, p1.z);
+               }
+               glEnd();
+       }
 }
 
 static void
@@ -982,7 +1039,7 @@ temporarySelection(MainView *mview, int flags, int clickCount, int ignoreExpande
                                continue;
                        if (mview->draggingMode == kMainViewSelectingRegion) {
                                /*  Check if this atom is within the selection rectangle  */
-                               GLfloat objectPos[3];
+                               double objectPos[3];
                                GLfloat screenPos[3];
                                Vector r1;
                                r1 = ap->r;
@@ -1010,80 +1067,85 @@ drawGraphite(MainView *mview)
        static GLfloat sDarkCyanColor[] = {0, 0.75, 0.75, 1};
        MDArena *arena;
        MDGraphiteArena *graphite;
+       Vector xaxis, yaxis, zaxis, origin;
+       Double R;
+       int i, j, i0, i1, j0, j1, ir;
+       Double x, dx, y, dy, xx, yy;
+       GLfloat p[12];
        if ((arena = mview->mol->arena) != NULL && (graphite = arena->graphite) != NULL) {
-               Vector xaxis, yaxis, zaxis, origin;
-               Double R;
-               int i, j, i0, i1, j0, j1, ir;
-               Double x, dx, y, dy, xx, yy;
-               GLfloat p[12];
                graphite_get_axes(graphite, &origin, &xaxis, &yaxis, &zaxis, &R);
-               i0 = -(mview->showGraphite / 2) - 1;
-               i1 = i0 + mview->showGraphite + 1;
-               j0 = -(mview->showGraphite / 2);
-               j1 = j0 + mview->showGraphite;
-               dx = 0.5 * R;
-               dy = 0.866025403784439 * R;
-               glDisable(GL_LIGHTING); 
-               glColor3fv(sDarkCyanColor);
-               for (i = i0; i <= i1; i++) {
-                       for (j = j0; j <= j1; j++) {
-                               Byte f1, f2, f3;
-                               ir = (i % 2 == 0 ? 0 : 1);
-                               x = 3 * i * dx;
-                               y = (2 * j + ir) * dy;
-                               yy = y - dy;
-                               xx = x - 2 * dx;
-                               p[0] = xaxis.x * xx + yaxis.x * y + origin.x;
-                               p[1] = xaxis.y * xx + yaxis.y * y + origin.y;
-                               p[2] = xaxis.z * xx + yaxis.z * y + origin.z;
-                               xx += dx;
-                               p[3] = xaxis.x * xx + yaxis.x * yy + origin.x;
-                               p[4] = xaxis.y * xx + yaxis.y * yy + origin.y;
-                               p[5] = xaxis.z * xx + yaxis.z * yy + origin.z;
-                               xx += 2 * dx;
-                               p[6] = xaxis.x * xx + yaxis.x * yy + origin.x;
-                               p[7] = xaxis.y * xx + yaxis.y * yy + origin.y;
-                               p[8] = xaxis.z * xx + yaxis.z * yy + origin.z;
-                               xx += dx;
-                               p[9] = xaxis.x * xx + yaxis.x * y + origin.x;
-                               p[10] = xaxis.y * xx + yaxis.y * y + origin.y;
-                               p[11] = xaxis.z * xx + yaxis.z * y + origin.z;
-                               f1 = f2 = f3 = 1;
-                               if (i == i0) {
-                                       f1 = f2 = 0;
-                                       if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
-                                               continue;
-                               } else if (i == i1) {
-                                       f2 = f3 = 0;
-                                       if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
-                                               continue;
-                               } else if (j == j1) {
-                                       if (ir == 1) {
-                                               f1 = f3 = 0;
-                                       } else if (i == i0 + 1) {
-                                               f1 = 0;
-                                       } else if (i == i1 - 1) {
-                                               f3 = 0;
-                                       }
-                               }
-                               glBegin(GL_LINES);
-                               if (f1) { 
-                                       glVertex3fv(p);
-                                       glVertex3fv(p + 3);
-                               }
-                               if (f2) {
-                                       glVertex3fv(p + 3);
-                                       glVertex3fv(p + 6);
-                               }
-                               if (f3) {
-                                       glVertex3fv(p + 6);
-                                       glVertex3fv(p + 9);
+       } else {
+               origin.x = origin.y = origin.z = 0.0;
+               xaxis.x = yaxis.y = zaxis.z = 1.0;
+               xaxis.y = xaxis.z = yaxis.x = yaxis.z = zaxis.x = zaxis.y = 0.0;
+               R = 1.42;
+       }
+       i0 = -(mview->showGraphite / 2) - 1;
+       i1 = i0 + mview->showGraphite + 1;
+       j0 = -(mview->showGraphite / 2);
+       j1 = j0 + mview->showGraphite;
+       dx = 0.5 * R;
+       dy = 0.866025403784439 * R;
+       glDisable(GL_LIGHTING); 
+       glColor3fv(sDarkCyanColor);
+       for (i = i0; i <= i1; i++) {
+               for (j = j0; j <= j1; j++) {
+                       Byte f1, f2, f3;
+                       ir = (i % 2 == 0 ? 0 : 1);
+                       x = 3 * i * dx;
+                       y = (2 * j + ir) * dy;
+                       yy = y - dy;
+                       xx = x - 2 * dx;
+                       p[0] = xaxis.x * xx + yaxis.x * y + origin.x;
+                       p[1] = xaxis.y * xx + yaxis.y * y + origin.y;
+                       p[2] = xaxis.z * xx + yaxis.z * y + origin.z;
+                       xx += dx;
+                       p[3] = xaxis.x * xx + yaxis.x * yy + origin.x;
+                       p[4] = xaxis.y * xx + yaxis.y * yy + origin.y;
+                       p[5] = xaxis.z * xx + yaxis.z * yy + origin.z;
+                       xx += 2 * dx;
+                       p[6] = xaxis.x * xx + yaxis.x * yy + origin.x;
+                       p[7] = xaxis.y * xx + yaxis.y * yy + origin.y;
+                       p[8] = xaxis.z * xx + yaxis.z * yy + origin.z;
+                       xx += dx;
+                       p[9] = xaxis.x * xx + yaxis.x * y + origin.x;
+                       p[10] = xaxis.y * xx + yaxis.y * y + origin.y;
+                       p[11] = xaxis.z * xx + yaxis.z * y + origin.z;
+                       f1 = f2 = f3 = 1;
+                       if (i == i0) {
+                               f1 = f2 = 0;
+                               if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
+                                       continue;
+                       } else if (i == i1) {
+                               f2 = f3 = 0;
+                               if ((ir == 0 && j == j0) || (ir == 1 && j == j1))
+                                       continue;
+                       } else if (j == j1) {
+                               if (ir == 1) {
+                                       f1 = f3 = 0;
+                               } else if (i == i0 + 1) {
+                                       f1 = 0;
+                               } else if (i == i1 - 1) {
+                                       f3 = 0;
                                }
-                               glEnd();
                        }
+                       glBegin(GL_LINES);
+                       if (f1) { 
+                               glVertex3fv(p);
+                               glVertex3fv(p + 3);
+                       }
+                       if (f2) {
+                               glVertex3fv(p + 3);
+                               glVertex3fv(p + 6);
+                       }
+                       if (f3) {
+                               glVertex3fv(p + 6);
+                               glVertex3fv(p + 9);
+                       }
+                       glEnd();
                }
-               glEnable(GL_LIGHTING);  
        }
+       enableLighting();
 }
 
 static GLfloat sRedColor[] = {1, 0, 0, 1};
@@ -1097,6 +1159,8 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
        int expanded = 0;
        Vector r1;
        GLfloat p[6];
+       double pp[3];
+       double rad;
        char label[16];
        GLfloat rgba[4];
        Transform *trp = NULL;
@@ -1105,26 +1169,12 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
                /*  Extra 2 atoms for the bond being newly created  */
                if (mview->draggingMode != kMainViewCreatingBond)
                        return;
-       /*      printf("mview->tempAtoms[%d] = %d\n", i - natoms, mview->tempAtoms[i - natoms]); */
                if (mview->tempAtoms[i1 - natoms] >= 0)
                        return;  /*  Already drawn  */
                ap = NULL;
                an1 = 6;
                r1 = mview->tempAtomPos[i1 - natoms];
                label[0] = 0;
-       } else if (i1 < 0) {
-               ExAtom *ep = mview->mol->exatoms + (- i1 - 1);
-               ap = ATOM_AT_INDEX(mview->mol->atoms, ep->index);
-               an1 = ap->atomicNumber;
-               r1 = ap->r;
-               trp = &(mview->mol->syms[ep->symop]);
-               if (/* !mview->mol->is_xtal_coord && */ mview->mol->cell != NULL) {
-                       TransformVec(&r1, mview->mol->cell->rtr, &r1);
-                       TransformVec(&r1, *trp, &r1);
-                       TransformVec(&r1, mview->mol->cell->tr, &r1);
-               } else TransformVec(&r1, *trp, &r1);
-               VecInc(r1, ep->dr);
-               label[0] = 0;
        } else {
                ap = ATOM_AT_INDEX(mview->mol->atoms, i1);
                if (ap == NULL)
@@ -1150,7 +1200,10 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
        if (selected) {
                memcpy(rgba, sRedColor, sizeof(rgba));
        } else {
-               rgba[0] = dp->r; rgba[1] = dp->g; rgba[2] = dp->b; rgba[3] = 1.0;
+               rgba[0] = dp->red / 65535.0;
+               rgba[1] = dp->green / 65535.0;
+               rgba[2] = dp->blue / 65535.0;
+               rgba[3] = 1.0;
        }
        if (expanded || periodicOffset != NULL) {
                rgba[0] *= 0.5;
@@ -1160,18 +1213,15 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
        glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
        if (periodicOffset != NULL)
                VecInc(r1, *periodicOffset);
-/*     if (mview->mol->is_xtal_coord)
-               TransformVec(&r1, mview->mol->cell->tr, &r1); */
        p[0] = r1.x; p[1] = r1.y; p[2] = r1.z;
        if (mview->draggingMode == kMainViewDraggingSelectedAtoms && selected) {
                p[0] += dragOffset->x;
                p[1] += dragOffset->y;
                p[2] += dragOffset->z;
        }
+       rad = VDW_RADIUS(dp) * mview->atomRadius;
        if (mview->showEllipsoids) {
                if (ap != NULL && ap->aniso != NULL) {
-               /*      Double xp[3][3]; */
-               /*      int i, j; */
                        GLfloat elip[9];
                        Mat33 pmat2;
                        int i;
@@ -1179,39 +1229,28 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
                                MatrixMul(pmat2, mview->mol->cell->rtr, ap->aniso->pmat);
                                MatrixMul(pmat2, *((Mat33 *)trp), pmat2);
                                MatrixMul(pmat2, mview->mol->cell->tr, pmat2);
-                               MatrixTranspose(pmat2, pmat2);
+                               for (i = 0; i < 9; i++)
+                                       elip[i] = pmat2[i] * mview->probabilityScale;
                        } else {
-                               MatrixTranspose(pmat2, ap->aniso->pmat);
+                               for (i = 0; i < 9; i++)
+                                       elip[i] = ap->aniso->pmat[i] * mview->probabilityScale;
                        }
-                       for (i = 0; i < 9; i++)
-                               elip[i] = pmat2[i] * mview->probabilityScale;
-               /*      for (i = 0; i < 3; i++) {
-                               Double w = ap->aniso->val[i];
-                               Vector vv = ap->aniso->axis[i];
-                               if (w <= 0.0)
-                                       w = 0.001;
-                               if (trp != NULL) {
-                                       TransformVec(&vv, mview->mol->cell->rtr, &vv);
-                                       MatrixVec(&vv, *((Mat33 *)trp), &vv);
-                                       TransformVec(&vv, mview->mol->cell->tr, &vv);
-                               }
-                               w *= mview->probabilityScale;
-                               xp[i][0] = w * vv.x;
-                               xp[i][1] = w * vv.y;
-                               xp[i][2] = w * vv.z;
-                       } */
-                       
-                       drawEllipsoid(p, elip, elip+3, elip+6, 15);
+                       drawEllipsoid(p, elip, elip+3, elip+6, mview->atomResolution * 3 / 2); /* Use higher resolution than spheres */
                } else {
-                       Double rad;
-                       rad = biso2radius(ap->tempFactor);
-                       rad *= mview->probabilityScale;
-                       drawSphere(p, rad, 8);
+                       if (ap != NULL) {
+                               //  Recalculate radius from temperature factor
+                               rad = biso2radius(ap->tempFactor);
+                               rad *= mview->probabilityScale;
+                       }
+                       drawSphere(p, rad, mview->atomResolution);
                }
        } else {
-               drawSphere(p, dp->radius * mview->atomRadius, 8);
+               drawSphere(p, rad, mview->atomResolution);
        }
-       if (MainView_convertObjectPositionToScreenPosition(mview, p, p + 3)) {
+       pp[0] = p[0];
+       pp[1] = p[1];
+       pp[2] = p[2];
+       if (MainView_convertObjectPositionToScreenPosition(mview, pp, p + 3)) {
        /*      fprintf(stderr, "atom %d: {%f, %f, %f}\n", i1, p[3], p[4], p[5]); */
                float fp[3];
                fp[0] = p[3]; fp[1] = p[4]; fp[2] = p[5];
@@ -1220,17 +1259,21 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
 }
 
 static void
-drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft, const Vector *dragOffset, const Vector *periodicOffset)
+drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft, const Vector *dragOffset, const Vector *periodicOffset, int isAnchorBond)
 {
        const ElementPar *dp;
        int i, in;
        int an[2];
-       int expanded[2];
+       char expanded[2];
+       char anchor[2];
        Vector r[2];
        GLfloat p[6];
        GLfloat rgba[4];
+       float rad_mul = 1.0;
+       float alpha_mul = 1.0;
        int natoms = mview->mol->natoms;
        expanded[0] = expanded[1] = 0;
+       anchor[0] = anchor[1] = 0;
 
        for (i = 0; i < 2; i++) {
                const Atom *ap;
@@ -1245,21 +1288,14 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
                                r[i] = mview->tempAtomPos[in - natoms];
                                an[i] = 6;
                        }
-               } else if (in < 0) {
-                       ExAtom *ep = mview->mol->exatoms + (- in - 1);
-                       if (!mview->showExpandedAtoms)
-                               return;
-                       ap = ATOM_AT_INDEX(mview->mol->atoms, ep->index);
-                       an[i] = ap->atomicNumber;
-                       r[i] = ap->r;
-                       TransformVec(&r[i], mview->mol->syms[ep->symop], &r[i]);
-                       VecInc(r[i], ep->dr);
                } else {
                        ap = ATOM_AT_INDEX(mview->mol->atoms, in);
                        an[i] = ap->atomicNumber;
                        r[i] = ap->r;
                        if (SYMOP_ALIVE(ap->symop))
                                expanded[i] = 1;
+                       if (ap->anchor != NULL)
+                               anchor[i] = 1;
                }
                if (!mview->showHydrogens && an[i] == 1)
                        return;
@@ -1275,19 +1311,24 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
                VecInc(r[0], *periodicOffset);
                VecInc(r[1], *periodicOffset);
        }
-/*     if (mview->mol->is_xtal_coord) {
-               TransformVec(&r[0], mview->mol->cell->tr, &r[0]);
-               TransformVec(&r[1], mview->mol->cell->tr, &r[1]);
-       } */
 
+       if (anchor[0] + anchor[1] == 2)
+               alpha_mul = 0.5;
+       if (anchor[0] + anchor[1] == 1)
+               rad_mul = (isAnchorBond ? 0.3 : 0.6);
+       
        dp = &(gElementParameters[an[0]]);
        if (dp == NULL)
                return;
        if (selected && selected2) {
                memcpy(rgba, sRedColor, sizeof(rgba));
        } else {
-               rgba[0] = dp->r; rgba[1] = dp->g; rgba[2] = dp->b; rgba[3] = 1.0;
+               rgba[0] = dp->red / 65535.0;
+               rgba[1] = dp->green / 65535.0;
+               rgba[2] = dp->blue / 65535.0;
+               rgba[3] = 1.0;
        }
+       rgba[3] *= alpha_mul;
        if (expanded[0] || periodicOffset != NULL) {
                rgba[0] *= 0.5;
                rgba[1] *= 0.5;
@@ -1311,14 +1352,18 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
                glEnd();
        } else {
                glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
-               drawCylinder(p, p + 3, mview->bondRadius, 6);
+               drawCylinder(p, p + 3, mview->bondRadius * rad_mul, mview->bondResolution, 0);
        }
        dp = &(gElementParameters[an[1]]);
        if (dp == NULL)
                return;
        if (!selected || !selected2) {
-               rgba[0] = dp->r; rgba[1] = dp->g; rgba[2] = dp->b; rgba[3] = 1.0;
+               rgba[0] = dp->red / 65535.0;
+               rgba[1] = dp->green / 65535.0;
+               rgba[2] = dp->blue / 65535.0;
+               rgba[3] = 1.0;
        }
+       rgba[3] *= alpha_mul;
        if (expanded[1] || periodicOffset != NULL) {
                rgba[0] *= 0.5;
                rgba[1] *= 0.5;
@@ -1334,7 +1379,7 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
                glEnd();
        } else {
                glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
-               drawCylinder(p, p + 3, mview->bondRadius, 6);
+               drawCylinder(p, p + 3, mview->bondRadius * rad_mul, 8, 0);
        }
 }
 
@@ -1342,7 +1387,7 @@ drawBond(MainView *mview, int i1, int i2, int selected, int selected2, int draft
 static void
 calcDragOffset(MainView *mview, Vector *outVector)
 {
-       GLfloat p[6];
+       double p[6];
        if (mview->draggingMode == kMainViewDraggingSelectedAtoms
        && MainView_convertScreenPositionToObjectPosition(mview, mview->dragStartPos, p)
        && MainView_convertScreenPositionToObjectPosition(mview, mview->dragEndPos, p + 3)) {
@@ -1356,6 +1401,32 @@ calcDragOffset(MainView *mview, Vector *outVector)
 }
 
 static void
+drawSurface(MainView *mview)
+{
+       int i, sn, k;
+       GLfloat rgba[4];
+       MCube *mc;
+       if (mview->mol == NULL || mview->mol->mcube == NULL || mview->mol->mcube->hidden != 0)
+               return;
+       mc = mview->mol->mcube;
+       for (sn = 0; sn <= 1; sn++) {
+               if (mc->c[sn].ntriangles == 0)
+                       continue;
+               for (i = 0; i < 4; i++)
+                       rgba[i] = mc->c[sn].rgba[i];
+               k = (sn == 0 ? -1 : 1);
+               glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, rgba);
+               glBegin(GL_TRIANGLES);
+               for (i = 0; mc->c[sn].triangles[i] >= 0; i++) {
+                       MCubePoint *mcp = mc->c[sn].cubepoints + mc->c[sn].triangles[i];
+                       glNormal3f(mcp->grad[0] * k, mcp->grad[1] * k, mcp->grad[2] * k);
+                       glVertex3f(mcp->pos[0], mcp->pos[1], mcp->pos[2]);
+               }
+               glEnd();                
+       }
+}
+
+static void
 drawModel(MainView *mview)
 {
        Molecule *mol;
@@ -1363,11 +1434,11 @@ drawModel(MainView *mview)
        int amin, amax, bmin, bmax, cmin, cmax, da, db, dc;
        Byte original;
        Double atomRadius, bondRadius;
-       Vector dragOffset;
        Vector periodicOffset;
        Vector *axes;
        int selected, selected2;
        char *selectFlags;
+       Atom *ap;
        int draft = mview->lineMode;
 /*    static Double gray[] = {0.8, 0.8, 0.8, 1}; */
        
@@ -1380,28 +1451,28 @@ drawModel(MainView *mview)
     }
 */
        if (mview->draggingMode == kMainViewSelectingRegion)
-               selectFlags = temporarySelection(mview, mview->modifierFlags, 0, 1);
+               selectFlags = temporarySelection(mview, mview->modifierFlags, 0, 0);
        else selectFlags = NULL;
        
        if (mview->draggingMode == kMainViewDraggingSelectedAtoms)
-               calcDragOffset(mview, &dragOffset);
-       else dragOffset.x = dragOffset.y = dragOffset.z = 0;
+               calcDragOffset(mview, &mview->dragOffset);
+       else mview->dragOffset.x = mview->dragOffset.y = mview->dragOffset.z = 0;
        
-       if (mview->showGraphite) {
+       if (mview->showGraphite > 0 && mview->showGraphiteFlag) {
                drawGraphite(mview);
        }
        
        amin = amax = bmin = bmax = cmin = cmax = 0;
        if (mview->showExpandedAtoms && mol->cell != NULL) {
-               if (mol->cell->flags[0]) {
+               if (mol->cell->flags[0] && mview->showPeriodicImageFlag) {
                        amin = mview->showPeriodicImage[0];
                        amax = mview->showPeriodicImage[1];
                }
-               if (mol->cell->flags[1]) {
+               if (mol->cell->flags[1] && mview->showPeriodicImageFlag) {
                        bmin = mview->showPeriodicImage[2];
                        bmax = mview->showPeriodicImage[3];
                }
-               if (mol->cell->flags[2]) {
+               if (mol->cell->flags[2] && mview->showPeriodicImageFlag) {
                        cmin = mview->showPeriodicImage[4];
                        cmax = mview->showPeriodicImage[5];
                }
@@ -1420,7 +1491,7 @@ drawModel(MainView *mview)
                                                VecScaleInc(periodicOffset, axes[1], db);
                                                VecScaleInc(periodicOffset, axes[2], dc);
                                        }
-                                       for (i = 0; i < natoms; i++) {
+                                       for (i = 0, ap = mview->mol->atoms; i < natoms; i++, ap = ATOM_NEXT(ap)) {
                                                if (mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
                                                        /*  Mouse event is detected  */
                                                        draft = 1;
@@ -1432,28 +1503,41 @@ drawModel(MainView *mview)
                                                        selected = selectFlags[i];
                                                else
                                                        selected = MoleculeIsAtomSelected(mview->mol, i);
-                                               drawAtom(mview, i, selected, &dragOffset, (original ? NULL : &periodicOffset));
+                                               drawAtom(mview, i, selected, &mview->dragOffset, (original ? NULL : &periodicOffset));
+                                               if (ap->anchor != NULL) {
+                                                       /*  Draw anchor bonds  */
+                                                       Int j, k;
+                                                       Int *cp = AtomConnectData(&ap->anchor->connect);
+                                                       for (j = 0; j < ap->anchor->connect.count; j++) {
+                                                               k = cp[j];
+                                                               if (selectFlags != NULL)
+                                                                       selected2 = selectFlags[k];
+                                                               else
+                                                                       selected2 = MoleculeIsAtomSelected(mview->mol, k);
+                                                               drawBond(mview, i, k, selected, selected2, draft, &mview->dragOffset, (original ? NULL : &periodicOffset), 1);
+                                                       }
+                                               }
                                        }
                                }
        
                                if (draft == 0) {
                                        if (original) {
                                                /*  Extra atoms  */
-                                               drawAtom(mview, natoms, 1, &dragOffset, NULL);
-                                               drawAtom(mview, natoms + 1, 1, &dragOffset, NULL);
+                                               drawAtom(mview, natoms, 1, &mview->dragOffset, NULL);
+                                               drawAtom(mview, natoms + 1, 1, &mview->dragOffset, NULL);
                                        }
                                        /*  Expanded atoms  */
-                                       if (mview->showExpandedAtoms) {
+                               /*      if (mview->showExpandedAtoms) {
                                                for (i = 0; i < mview->mol->nexatoms; i++) {
                                                        if (mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
-                                                               /*  Mouse event is detected  */
+                                                               //  Mouse event is detected
                                                                draft = 1;
                                                                break;
-                                                       /*      goto cleanup;  */
+                                                       //      goto cleanup;
                                                        }
                                                        drawAtom(mview, -i-1, (selectFlags == NULL ? 0 : selectFlags[i]), &dragOffset, (original ? NULL : &periodicOffset));
                                                }
-                                       }
+                                       } */
                                }
                        }
                }
@@ -1490,22 +1574,22 @@ skip:
                                                selected = selectFlags[n1];
                                                selected2 = selectFlags[n2];
                                        }
-                                       drawBond(mview, n1, n2, selected, selected2, draft, &dragOffset, (original ? NULL : &periodicOffset));
+                                       drawBond(mview, n1, n2, selected, selected2, draft, &mview->dragOffset, (original ? NULL : &periodicOffset), 0);
                                }
                                
                                /*  Extra bond  */
                                if (original && mview->draggingMode == kMainViewCreatingBond) {
-                                       drawBond(mview, natoms, natoms + 1, 1, 1, draft, &dragOffset, NULL);
+                                       drawBond(mview, natoms, natoms + 1, 1, 1, draft, &mview->dragOffset, NULL, 0);
                                }
                                
                                /*  Expanded bonds  */
-                               for (i = 0; i < mview->mol->nexbonds; i++) {
+                       /*      for (i = 0; i < mview->mol->nexbonds; i++) {
                                        int n1, n2;
                                        if (draft == 0 && mview->draggingMode != 0 && i % 50 == 0 && MainViewCallback_mouseCheck(mview)) {
-                                               /*  Mouse event is detected  */
+                                               //  Mouse event is detected
                                                draft = 1;
                                                glDisable(GL_LIGHTING);
-                                       /*      goto cleanup;  */
+                                       //      goto cleanup;
                                        }
                                        n1 = mview->mol->exbonds[i * 2];
                                        n2 = mview->mol->exbonds[i * 2 + 1];
@@ -1521,118 +1605,194 @@ skip:
                                                selected2 = selectFlags[n2];
                                        }
                                        drawBond(mview, mview->mol->exbonds[i * 2], mview->mol->exbonds[i * 2 + 1], selected, selected2, draft, &dragOffset, (original ? NULL : &periodicOffset));
-                               }
+                               } */
                        }
                }
        }
        
 /*  cleanup: */
        if (draft)
-               glEnable(GL_LIGHTING);
+               enableLighting();
        if (selectFlags != NULL)
                free(selectFlags);
 }
 
 static void
+drawGraphics(MainView *mview)
+{
+       int i, j;
+       MainViewGraphic *g;
+       for (i = 0; i < mview->ngraphics; i++) {
+               g = &mview->graphics[i];
+               if (g->visible == 0)
+                       continue;
+               glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, g->rgba);
+               switch (g->kind) {
+                       case kMainViewGraphicLine:
+                               glDisable(GL_LIGHTING);
+                               glColor4fv(g->rgba);
+                               glBegin(GL_LINE_STRIP);
+                               for (j = 0; j < g->npoints; j++) {
+                                       if (g->points[j * 3] >= kInvalidFloat)
+                                               break;
+                                       glVertex3fv(&g->points[j * 3]);
+                               }
+                               glEnd();
+                               enableLighting();
+                               break;
+                       case kMainViewGraphicPoly: {
+                       /*      Vector v0, v1, v2, v3; */
+                               glBegin(GL_TRIANGLE_FAN);
+                       /*      v1.x = g->points[0] - g->points[g->npoints - 3];
+                               v1.y = g->points[1] - g->points[g->npoints - 2];
+                               v1.z = g->points[2] - g->points[g->npoints - 1];
+                               v0 = v1;  */
+                               for (j = 0; j < g->npoints; j++) {
+                               /*      v2.x = g->points[j * 3 + 3] - g->points[j * 3];
+                                       v2.y = g->points[j * 3 + 4] - g->points[j * 3 + 1];
+                                       v2.z = g->points[j * 3 + 5] - g->points[j * 3 + 2];
+                                       VecCross(v3, v1, v2);
+                                       if (NormalizeVec(&v3, &v3) == 0)
+                                               glNormal3f(v3.x, v3.y, v3.z); */
+                                       glNormal3fv(&g->normals[j * 3]);
+                                       glVertex3fv(&g->points[j * 3]);
+                               /*      v1 = v2; */
+                               }
+                               if (g->closed) {
+                               /*      VecCross(v3, v1, v0);
+                                       if (NormalizeVec(&v3, &v3) == 0)
+                                               glNormal3f(v3.x, v3.y, v3.z); */
+                                       glNormal3fv(g->normals);
+                                       glVertex3fv(g->points);
+                               }
+                               glEnd();
+                               break;
+                       }
+                       case kMainViewGraphicCylinder:
+                               drawCylinder(g->points, g->points + 3, g->points[6], 15, g->closed);
+                               break;
+                       case kMainViewGraphicCone:
+                               drawCone(g->points, g->points + 3, g->points[6], 15, g->closed);
+                               break;
+                       case kMainViewGraphicEllipsoid:
+                               drawEllipsoid(g->points, g->points + 3, g->points + 6, g->points + 9, 8);
+                               break;
+               }
+       }
+}
+
+static void
 drawUnitCell(MainView *mview)
 {
        GLfloat a[3], b[3], c[3], ab[3], bc[3], ca[3], abc[3];
        XtalCell *cp;
-       static GLfloat sGrayColor[] = {0.5, 0.5, 0.5, 1};
        GLfloat origin[3];
-       int i;
-       glDisable(GL_LIGHTING);
-       for (i = 0; i < 1; i++) {
-               if (i == 0) {
-                       if (!mview->showUnitCell || (cp = mview->mol->cell) == NULL)
-                               continue;
-                       a[0] = cp->tr[0];
-                       a[1] = cp->tr[3];
-                       a[2] = cp->tr[6];
-                       b[0] = cp->tr[1];
-                       b[1] = cp->tr[4];
-                       b[2] = cp->tr[7];
-                       c[0] = cp->tr[2];
-                       c[1] = cp->tr[5];
-                       c[2] = cp->tr[8];
-                       origin[0] = cp->origin.x;
-                       origin[1] = cp->origin.y;
-                       origin[2] = cp->origin.z;
-                       glColor3f(0.75, 0.2, 0.0);
-/*             } else {
-                       if (!mview->showPeriodicBox || (bp = mview->mol->box) == NULL)
-                               continue;
-                       origin[0] = bp->origin.x;
-                       origin[1] = bp->origin.y;
-                       origin[2] = bp->origin.z;
-                       a[0] = bp->axes[0].x + origin[0];
-                       a[1] = bp->axes[0].y + origin[1];
-                       a[2] = bp->axes[0].z + origin[2];
-                       b[0] = bp->axes[1].x + origin[0];
-                       b[1] = bp->axes[1].y + origin[1];
-                       b[2] = bp->axes[1].z + origin[2];
-                       c[0] = bp->axes[2].x + origin[0];
-                       c[1] = bp->axes[2].y + origin[1];
-                       c[2] = bp->axes[2].z + origin[2];
-                       glColor3f(0.0, 0.75, 0.2); */
-               }
-       
-               ab[0] = a[0] + b[0]; ab[1] = a[1] + b[1]; ab[2] = a[2] + b[2];
-               bc[0] = b[0] + c[0]; bc[1] = b[1] + c[1]; bc[2] = b[2] + c[2];
-               ca[0] = c[0] + a[0]; ca[1] = c[1] + a[1]; ca[2] = c[2] + a[2];
-               abc[0] = a[0] + bc[0]; abc[1] = a[1] + bc[1]; abc[2] = a[2] + bc[2];
-
-       /*      glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, sGrayColor);
-               drawCylinder(sOrigin, a, 0.05, 4);
-               drawCylinder(sOrigin, b, 0.05, 4);
-               drawCylinder(sOrigin, c, 0.05, 4);
-               drawCylinder(a, ab, 0.05, 4);
-               drawCylinder(a, ca, 0.05, 4);
-               drawCylinder(b, ab, 0.05, 4);
-               drawCylinder(b, bc, 0.05, 4);
-               drawCylinder(c, ca, 0.05, 4);
-               drawCylinder(c, bc, 0.05, 4);
-               drawCylinder(ab, abc, 0.05, 4);
-               drawCylinder(bc, abc, 0.05, 4);
-               drawCylinder(ca, abc, 0.05, 4); */
+       if (!mview->showUnitCell || (cp = mview->mol->cell) == NULL)
+               return;
+       origin[0] = cp->origin.x;
+       origin[1] = cp->origin.y;
+       origin[2] = cp->origin.z;
+       a[0] = cp->axes[0].x + origin[0];
+       a[1] = cp->axes[0].y + origin[1];
+       a[2] = cp->axes[0].z + origin[2];
+       b[0] = cp->axes[1].x + origin[0];
+       b[1] = cp->axes[1].y + origin[1];
+       b[2] = cp->axes[1].z + origin[2];
+       c[0] = cp->axes[2].x + origin[0];
+       c[1] = cp->axes[2].y + origin[1];
+       c[2] = cp->axes[2].z + origin[2];
+
+       ab[0] = a[0] + b[0]; ab[1] = a[1] + b[1]; ab[2] = a[2] + b[2];
+       bc[0] = b[0] + c[0]; bc[1] = b[1] + c[1]; bc[2] = b[2] + c[2];
+       ca[0] = c[0] + a[0]; ca[1] = c[1] + a[1]; ca[2] = c[2] + a[2];
+       abc[0] = a[0] + bc[0]; abc[1] = a[1] + bc[1]; abc[2] = a[2] + bc[2];
 
-               glBegin(GL_LINES);
-               glVertex3fv(origin);
-               glVertex3fv(a);
-               glVertex3fv(origin);
-               glVertex3fv(b);
-               glVertex3fv(origin);
-               glVertex3fv(c);
-               glVertex3fv(a);
-               glVertex3fv(ab);
-               glVertex3fv(a);
-               glVertex3fv(ca);
-               glVertex3fv(b);
-               glVertex3fv(ab);
-               glVertex3fv(b);
-               glVertex3fv(bc);
-               glVertex3fv(c);
-               glVertex3fv(ca);
-               glVertex3fv(c);
-               glVertex3fv(bc);
-               glVertex3fv(ab);
-               glVertex3fv(abc);
-               glVertex3fv(bc);
-               glVertex3fv(abc);
-               glVertex3fv(ca);
-               glVertex3fv(abc);
-               glEnd();
-       }
-       glEnable(GL_LIGHTING);
+       glDisable(GL_LIGHTING);
+       glColor3f(0.75, 0.2, 0.0);
+       glBegin(GL_LINES);
+       glVertex3fv(origin);
+       glVertex3fv(a);
+       glVertex3fv(origin);
+       glVertex3fv(b);
+       glVertex3fv(origin);
+       glVertex3fv(c);
+       glVertex3fv(a);
+       glVertex3fv(ab);
+       glVertex3fv(a);
+       glVertex3fv(ca);
+       glVertex3fv(b);
+       glVertex3fv(ab);
+       glVertex3fv(b);
+       glVertex3fv(bc);
+       glVertex3fv(c);
+       glVertex3fv(ca);
+       glVertex3fv(c);
+       glVertex3fv(bc);
+       glVertex3fv(ab);
+       glVertex3fv(abc);
+       glVertex3fv(bc);
+       glVertex3fv(abc);
+       glVertex3fv(ca);
+       glVertex3fv(abc);
+       glEnd();
+
+       enableLighting();
+}
+
+/*  For debugging surface view  */
+static void
+drawCubeBoundary(MainView *mview)
+{
+    MCube *mc;
+    GLfloat ox, oy, oz;
+    GLfloat px, py, pz;
+    if ((mc = mview->mol->mcube) == NULL || mc->showbox == 0)
+        return;
+    ox = mc->origin.x;
+    oy = mc->origin.y;
+    oz = mc->origin.z;
+    px = ox + mc->dx * mc->nx;
+    py = oy + mc->dy * mc->ny;
+    pz = oz + mc->dz * mc->nz;
+    glDisable(GL_LIGHTING);
+    glColor3f(0.0, 0.5, 0.75);
+    glBegin(GL_LINES);
+    glVertex3f(ox, oy, oz);
+    glVertex3f(px, oy, oz);
+    glVertex3f(ox, oy, oz);
+    glVertex3f(ox, py, oz);
+    glVertex3f(ox, oy, oz);
+    glVertex3f(ox, oy, pz);
+    glVertex3f(px, oy, oz);
+    glVertex3f(px, py, oz);
+    glVertex3f(px, oy, oz);
+    glVertex3f(px, oy, pz);
+    glVertex3f(ox, py, oz);
+    glVertex3f(px, py, oz);
+    glVertex3f(ox, py, oz);
+    glVertex3f(ox, py, pz);
+    glVertex3f(ox, oy, pz);
+    glVertex3f(px, oy, pz);
+    glVertex3f(ox, oy, pz);
+    glVertex3f(ox, py, pz);
+    glVertex3f(px, py, oz);
+    glVertex3f(px, py, pz);
+    glVertex3f(px, oy, pz);
+    glVertex3f(px, py, pz);
+    glVertex3f(ox, py, pz);
+    glVertex3f(px, py, pz);
+    glEnd();
+    
+    enableLighting();
 }
 
 static void
 drawRotationCenter(MainView *mview)
 {
        GLfloat ps[3], pe[3], col[4];
-       float fd[2];  /*  Fovy and distance  */
-       float tr[3];  /*  Translation  */
-       float r, rr;
+       double fd[2];  /*  Fovy and distance  */
+       double tr[3];  /*  Translation  */
+       double r, rr;
        if (mview == NULL || !mview->showRotationCenter)
                return;
        TrackballGetTranslate(mview->track, tr);
@@ -1656,7 +1816,7 @@ drawRotationCenter(MainView *mview)
        glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
        drawSphere(ps, rr, 8);
        drawSphere(pe, rr, 8);
-       drawCylinder(ps, pe, rr, 8);
+       drawCylinder(ps, pe, rr, 8, 0);
        ps[0] = tr[0];
        ps[1] = tr[1] - r;
        pe[0] = tr[0];
@@ -1665,7 +1825,7 @@ drawRotationCenter(MainView *mview)
        glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
        drawSphere(ps, rr, 8);
        drawSphere(pe, rr, 8);
-       drawCylinder(ps, pe, rr, 8);
+       drawCylinder(ps, pe, rr, 8, 0);
        ps[1] = tr[1];
        ps[2] = tr[2] - r;
        pe[1] = tr[1];
@@ -1674,7 +1834,7 @@ drawRotationCenter(MainView *mview)
        glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, col);
        drawSphere(ps, rr, 8);
        drawSphere(pe, rr, 8);
-       drawCylinder(ps, pe, rr, 8);
+       drawCylinder(ps, pe, rr, 8, 0);
 }
 
 static int
@@ -1691,7 +1851,7 @@ compareLabelByDepth(const void *ap, const void *bp)
 static void
 drawLabels(MainView *mview)
 {
-       Transform *trp;
+/*     Transform *trp; */
        Atom *ap;
        LabelRecord *lp;
        int i, nlabels;
@@ -1774,8 +1934,8 @@ drawLabels(MainView *mview)
 void
 MainView_drawModel(MainView *mview)
 {
-    float w[4], dimension, distance;
-       float frame[4], width, height;
+    double w[4], dimension, distance;
+       float frame[4], width, height, scale;
        GLdouble *pp;
        Transform mtr;
        
@@ -1785,27 +1945,41 @@ MainView_drawModel(MainView *mview)
                return;
 
     if (!mview->isInitialized) {
-        MainView_initializeOpenGLView(mview);
+        MainView_initializeOpenGL();
+        mview->view_scale = MainViewCallback_getContentScaleFactor(mview);
+               mview->isInitialized = 1;
     }
     
-    /*  Clear the buffer  */
-    glClearColor (0, 0, 0, 0);
-    glClear(GL_COLOR_BUFFER_BIT |
-            GL_DEPTH_BUFFER_BIT);
-
-       if (mview->mol == NULL)
-               return;
-
        dimension = mview->dimension;
-/*    dimension = [model dimension];  */
+
+       if (mview->offline_scale == 0.0)
+        scale = mview->view_scale;
+       else
+               scale = mview->offline_scale;
 
        MainViewCallback_frame(mview, frame);
        width = frame[2] - frame[0];
        height = frame[3] - frame[1];
 
+       if (mview->offline_width > 0)
+               width = mview->offline_width;
+       if (mview->offline_height > 0)
+               height = mview->offline_height;
+       
+       width *= scale;
+       height *= scale;
+
     glViewport(0, 0, width, height);
     
-    /*  Set up the projection  */
+       /*  Clear the buffer  */
+    glClearColor(mview->background_color[0], mview->background_color[1], mview->background_color[2], mview->background_color[3]);
+    glClear(GL_COLOR_BUFFER_BIT |
+            GL_DEPTH_BUFFER_BIT);
+       
+       if (mview->mol == NULL)
+               return;
+
+       /*  Set up the projection  */
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();
        TrackballGetPerspective(mview->track, w);
@@ -1814,7 +1988,7 @@ MainView_drawModel(MainView *mview)
     mview->perspective_vector[1] = width / height;
     mview->perspective_vector[2] = dimension;
     mview->perspective_vector[3] = distance + 200.0 * dimension;
-    gluPerspective(mview->perspective_vector[0], mview->perspective_vector[1], mview->perspective_vector[2], mview->perspective_vector[3]);
+    myGluPerspective(mview->perspective_vector[0], mview->perspective_vector[1], mview->perspective_vector[2], mview->perspective_vector[3]);
 
     /*  Set up the model view  */
     glMatrixMode(GL_MODELVIEW);
@@ -1833,21 +2007,24 @@ MainView_drawModel(MainView *mview)
        
        MainViewCallback_clearLabels(mview);
     drawModel(mview);
+       drawSurface(mview);
        drawUnitCell(mview);
+    drawCubeBoundary(mview);
        drawRotationCenter(mview);
-
+       drawGraphics(mview);
+       
        /*  Get important matrices and vectors  */
     glGetDoublev(GL_MODELVIEW_MATRIX, mview->modelview_matrix);
     glGetDoublev(GL_PROJECTION_MATRIX, mview->projection_matrix);
        pp = mview->modelview_matrix;
-       mtr[0] = pp[0]; mtr[1] = pp[4]; mtr[2] = pp[8];
-       mtr[3] = pp[1]; mtr[4] = pp[5]; mtr[5] = pp[9];
-       mtr[6] = pp[2]; mtr[7] = pp[6]; mtr[8] = pp[10];
+       mtr[0] = pp[0]; mtr[1] = pp[1]; mtr[2] = pp[2];
+       mtr[3] = pp[4]; mtr[4] = pp[5]; mtr[5] = pp[6];
+       mtr[6] = pp[8]; mtr[7] = pp[9]; mtr[8] = pp[10];
        mtr[9] = pp[12]; mtr[10] = pp[13]; mtr[11] = pp[14];
        TransformInvert(mtr, mtr);
        mview->camera.x = mtr[9]; mview->camera.y = mtr[10]; mview->camera.z = mtr[11];
-       mview->lookto.x = mtr[2]; mview->lookto.y = mtr[5]; mview->lookto.z = mtr[8];
-       mview->up.x = mtr[1]; mview->up.y = mtr[4]; mview->up.z = mtr[7];
+       mview->lookto.x = mtr[6]; mview->lookto.y = mtr[7]; mview->lookto.z = mtr[8];
+       mview->up.x = mtr[3]; mview->up.y = mtr[4]; mview->up.z = mtr[5];
 
        /*  Draw labels  */
        glDisable(GL_LIGHTING);
@@ -1861,14 +2038,14 @@ MainView_drawModel(MainView *mview)
        glLoadIdentity ();
        glMatrixMode (GL_MODELVIEW);
        glLoadIdentity ();
-       glOrtho(0, width, 0, -height, 0.0,-1.0);  /*  non-flipped view  */
+       glOrtho(0, width, 0, -height, 0.0, -1.0);  /*  non-flipped view  */
        drawLabels(mview);
 #if __WXMAC__
        glDisable (GL_TEXTURE_RECTANGLE_EXT);
 #endif
 //     glDisable(GL_BLEND);
 //     glEnable(GL_DEPTH_TEST);
-       glEnable(GL_LIGHTING);
+       enableLighting();
 
     if (mview->draggingMode == kMainViewSelectingRegion) {
                /*  Draw selection rectangle  */
@@ -1881,14 +2058,14 @@ MainView_drawModel(MainView *mview)
         glOrtho(0, width, 0, height, -1.0, 1.0);
         glColor3f(1.0, 1.0, 0.0);
         glBegin(GL_LINE_STRIP);
-               glVertex2f(mview->dragStartPos[0], mview->dragStartPos[1]);
-               glVertex2f(mview->dragStartPos[0], mview->dragEndPos[1]);
-               glVertex2f(mview->dragEndPos[0], mview->dragEndPos[1]);
-               glVertex2f(mview->dragEndPos[0], mview->dragStartPos[1]);
-               glVertex2f(mview->dragStartPos[0], mview->dragStartPos[1]);
+               glVertex2f(mview->dragStartPos[0] * scale, mview->dragStartPos[1] * scale);
+               glVertex2f(mview->dragStartPos[0] * scale, mview->dragEndPos[1] * scale);
+               glVertex2f(mview->dragEndPos[0] * scale, mview->dragEndPos[1] * scale);
+               glVertex2f(mview->dragEndPos[0] * scale, mview->dragStartPos[1] * scale);
+               glVertex2f(mview->dragStartPos[0] * scale, mview->dragStartPos[1] * scale);
         glEnd();
                glEnable(GL_DEPTH_TEST);
-               glEnable(GL_LIGHTING);
+               enableLighting();
     }
 
     glFinish();
@@ -2043,8 +2220,10 @@ showAtomsInInfoText(MainView *mview, int n1, int n2)
 void
 MainView_mouseDown(MainView *mview, const float *mousePos, int flags)
 {
-       float p[3];
-       float screenPos[3], objectPos[3];
+       double p[3];
+       float fp[3];
+       float screenPos[3];
+       double objectPos[3];
        int n1, n2, found;
        Atom *ap;
        mview->dragStartPos[0] = mousePos[0];
@@ -2072,11 +2251,11 @@ MainView_mouseDown(MainView *mview, const float *mousePos, int flags)
                                objectPos[0] = r2.x;
                                objectPos[1] = r2.y;
                                objectPos[2] = r2.z;
-                               if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, p)) {
+                               if (MainView_convertObjectPositionToScreenPosition(mview, objectPos, fp)) {
                                        Double w;
-                                       r1.x = p[0] - screenPos[0];
-                                       r1.y = p[1] - screenPos[1];
-                                       r1.z = p[2] - screenPos[2];
+                                       r1.x = fp[0] - screenPos[0];
+                                       r1.y = fp[1] - screenPos[1];
+                                       r1.z = fp[2] - screenPos[2];
                                        NormalizeVec(&r1, &r1);
                                        r2.x = mousePos[0] - screenPos[0];
                                        r2.y = mousePos[1] - screenPos[1];
@@ -2106,8 +2285,8 @@ MainView_mouseDown(MainView *mview, const float *mousePos, int flags)
                case kTrackballTranslateMode:
                case kTrackballScaleMode:
                        if (!found) {
-                               mousePosToTrackballPos(mview, mousePos, p);
-                               TrackballStartDragging(mview->track, p, mview->mode);
+                               mousePosToTrackballPos(mview, mousePos, fp);
+                               TrackballStartDragging(mview->track, fp, mview->mode);
                                mview->draggingMode = kMainViewMovingTrackball;
                                break;
                        } else {
@@ -2196,7 +2375,7 @@ MainView_mouseDragged(MainView *mview, const float *mousePos, int flags)
 {
        float p[2];
        if (mview->isDragging == 0) {
-               if (abs(mousePos[0] - mview->dragStartPos[0]) >= 3 || abs(mousePos[1] - mview->dragStartPos[1]) >= 3)
+               if (fabsf(mousePos[0] - mview->dragStartPos[0]) >= 3 || fabsf(mousePos[1] - mview->dragStartPos[1]) >= 3)
                        mview->isDragging = 1;
                else return;
        }
@@ -2222,7 +2401,8 @@ MainView_mouseDragged(MainView *mview, const float *mousePos, int flags)
                                mview->tempAtoms[1] = n1;
                                mview->tempAtomPos[1] = ap->r;
                        } else {
-                               float screenPos[3], objectPos[3];
+                               float screenPos[3];
+                               double objectPos[3];
                                Vector r1;
                                mview->tempAtoms[1] = -1;
                                /*  Convert the position of temporary atom 0 to screen position */
@@ -2311,10 +2491,20 @@ MainView_mouseUp(MainView *mview, const float *mousePos, int flags, int clickCou
 
        if (clickCount == 2 && mview->isDragging == 0) {
                /*  Create a new molecular fragment  */
-               buf[0] = 0;
+               int n1, n2, status, found;
+        char *p;
+               found = MainView_findObjectAtPoint(mview, mousePos, &n1, &n2, 0, 0);
+        status = MyAppCallback_getGlobalSettingsWithType("global.entered_formula", 's', &p);
+        if (status == 0) {
+            strncpy(buf, p, sizeof buf - 1);
+            buf[sizeof buf - 1] = 0;
+        } else {
+            buf[0] = 0;
+        }
                if (MyAppCallback_getTextWithPrompt("Enter formula (e.g. CH2OCH3)", buf, sizeof buf) == 0)
                        return;
-               MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("s"), "dock_formula", buf);
+        MyAppCallback_setGlobalSettingsWithType("global.entered_formula", 's', buf);
+               MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("si"), "cmd_fuse_or_dock", buf, n1);
                goto exit;
        }
 
@@ -2335,10 +2525,14 @@ MainView_mouseUp(MainView *mview, const float *mousePos, int flags, int clickCou
                                MolActionCreateAndPerform(mview->mol, gMolActionDeleteAnAtom, (Int)n1);
                                
                        } else {
-                               Int nn[2];
-                               nn[0] = n1;
-                               nn[1] = n2;
-                               MolActionCreateAndPerform(mview->mol, gMolActionDeleteBonds, 2, nn);
+                               Int bn = MoleculeLookupBond(mview->mol, n1, n2);
+                               if (bn >= 0) {
+                                       IntGroup *ig = IntGroupNewWithPoints(bn, 1, -1);
+                                       MolActionCreateAndPerform(mview->mol, gMolActionDeleteBonds, ig);
+                                       IntGroupRelease(ig);
+                               } else {
+                                       fprintf(stderr, "Internal error %s:%d: bond to delete is not found", __FILE__, __LINE__);
+                               }
                        }
                }
                goto exit;
@@ -2351,21 +2545,17 @@ MainView_mouseUp(MainView *mview, const float *mousePos, int flags, int clickCou
                case kMainViewDraggingSelectedAtoms: {
                        Vector offset;
                        if (mview->isDragging) {
+                               IntGroup *dupSelection = IntGroupNewFromIntGroup(mview->mol->selection);
                                calcDragOffset(mview, &offset);
                        /*      if (mview->mol->is_xtal_coord)
                                        TransformVec(&offset, mview->mol->cell->rtr, &offset); */
-                               MolActionCreateAndPerform(mview->mol, gMolActionTranslateAtoms, &offset, mview->mol->selection);
-               /*      } else if (clickCount == 2) {
-                               buf[0] = 0;
-                               if (MyAppCallback_getTextWithPrompt("Enter formula (e.g. CH2OCH3)", buf, sizeof buf) == 0)
-                                       return;
-                               MolActionCreateAndPerform(mview->mol, SCRIPT_ACTION("s"), "dock_formula", buf);
-               */
+                               MolActionCreateAndPerform(mview->mol, gMolActionTranslateAtoms, &offset, dupSelection);
+                               IntGroupRelease(dupSelection);
                        }
                        break;
                }
                case kMainViewSelectingRegion: {
-                       char *selectFlags = temporarySelection(mview, mview->modifierFlags, (mview->isDragging ? 0 : 1), 1);
+                       char *selectFlags = temporarySelection(mview, mview->modifierFlags, (mview->isDragging ? 0 : 1), 0);
                        if (selectFlags != NULL) {
                                IntGroup *ig = IntGroupNew();
                                if (ig != NULL) {
@@ -2416,7 +2606,7 @@ MainView_mouseUp(MainView *mview, const float *mousePos, int flags, int clickCou
                                b[0] = n1;
                                b[1] = n2;
                                b[2] = kInvalidIndex;
-                               MolActionCreateAndPerform(mview->mol, gMolActionAddBonds, 2, b);
+                               MolActionCreateAndPerform(mview->mol, gMolActionAddBonds, 2, b, NULL);
                        /*      MoleculeAddBonds(mview->mol, b, NULL); */
                        }
                        break;
@@ -2453,7 +2643,7 @@ MainView_rotateBySlider(MainView *mview, float angle, int mode, int mouseStatus,
                        } else return;
                } else if ((modifierFlags & kAltKeyMask) != 0) {
                        /*  Rotate selection  */
-                       float tr[3];
+                       double tr[3];
                        if (mview->mol->selection == NULL)
                                return;
                        mview->rotateFragment = IntGroupNewFromIntGroup(mview->mol->selection);
@@ -2511,9 +2701,9 @@ MainView_rotateBySlider(MainView *mview, float angle, int mode, int mouseStatus,
                        }
                        
                } else {
-                       float quat[4];
-                       float cs = cos(angle / 2);
-                       float sn = sin(angle / 2);
+                       double quat[4];
+                       double cs = cos(angle / 2);
+                       double sn = sin(angle / 2);
                        if (mouseStatus == 0) { /* mouseUp */
                                TrackballEndDragging(mview->track, NULL);
                        } else { /* mouseDragged */
@@ -2584,7 +2774,7 @@ MainView_centerSelection(MainView *mview)
 {
        IntGroup *ig;
        Vector c;
-       float tr[3];
+       double tr[3];
        if (mview == NULL || mview->mol == NULL)
                return;
        ig = MoleculeGetSelection(mview->mol);
@@ -2611,7 +2801,7 @@ MainView_copy(MainView *mview)
        || mol2 == NULL
        || (p = MoleculeSerialize(mol2, &len, &time)) == NULL)
                return 2;
-       result = MoleculeCallback_writeToPasteboard(kMoleculePasteboardType, p, len);
+       result = MoleculeCallback_writeToPasteboard(gMoleculePasteboardType, p, len);
        if (result != 0)
                return result;
        if (mol2 != NULL)
@@ -2653,7 +2843,7 @@ MainView_cut(MainView *mview)
 int
 MainView_isPastable(MainView *mview)
 {
-       if (MoleculeCallback_isDataInPasteboard(kMoleculePasteboardType))
+       if (MoleculeCallback_isDataInPasteboard(gMoleculePasteboardType))
                return 1;
        else return 0;
 }
@@ -2672,7 +2862,7 @@ MainView_paste(MainView *mview)
        }       
        if (!MainView_isPastable(mview))
                return 1;
-       if (MoleculeCallback_readFromPasteboard(kMoleculePasteboardType, &p, &len) != 0)
+       if (MoleculeCallback_readFromPasteboard(gMoleculePasteboardType, &p, &len) != 0)
                return 2;
        if ((mol2 = MoleculeDeserialize(p, len, &time)) == NULL) {
                free(p);
@@ -2704,7 +2894,7 @@ MainView_paste(MainView *mview)
                        }
                }
                if (IntGroupGetCount(sel) > 0)
-                       MoleculeUnmerge(mol2, NULL, sel, 0);
+                       MoleculeUnmerge(mol2, NULL, sel, 0, NULL, NULL, 0);
                IntGroupRelease(sel);
        }
        
@@ -2748,11 +2938,13 @@ MainView_pasteParameters(MainView *mview)
        char *p, *pp;
        Int len, i;
        IntGroup *newsel;
-       if (mview == NULL || mview->mol == NULL || mview->mol->par == NULL)
+       if (mview == NULL || mview->mol == NULL)
                return -1;
-       if (!MoleculeCallback_isDataInPasteboard(kParameterPasteboardType))
+       if (mview->mol->par == NULL)
+               mview->mol->par = ParameterNew();
+       if (!MoleculeCallback_isDataInPasteboard(gParameterPasteboardType))
                return 1;
-       if (MoleculeCallback_readFromPasteboard(kParameterPasteboardType, (void **)&p, &len) != 0)
+       if (MoleculeCallback_readFromPasteboard(gParameterPasteboardType, (void **)&p, &len) != 0)
                return 2;
        pp = p;
        newsel = IntGroupNew();
@@ -2810,7 +3002,9 @@ MainView_copyOrCutParameters(MainView *mview, int flags)
        Parameter *par = mview->mol->par;
        void *pb_ptr = NULL;
        Int pb_len = 0;
-       i = 0;
+       if (ig == NULL || (i = IntGroupGetCount(ig)) == 0)
+               return 0;
+       i--;
        ig2 = NULL;
        while (1) {
                n = IntGroupGetNthPoint(ig, i);
@@ -2852,14 +3046,14 @@ MainView_copyOrCutParameters(MainView *mview, int flags)
                                ig2 = IntGroupNew();
                        IntGroupAdd(ig2, idx, 1);
                }
-               i++;
+               i--;
        }
        if (ig2 != NULL)
                IntGroupRelease(ig2);
        IntGroupRelease(ig);
        
        if (flags & 2) {
-               n = MoleculeCallback_writeToPasteboard(kParameterPasteboardType, pb_ptr, pb_len);
+               n = MoleculeCallback_writeToPasteboard(gParameterPasteboardType, pb_ptr, pb_len);
                if (n != 0)
                        return n;
        }
@@ -2884,32 +3078,48 @@ typedef struct ColumnInfoRecord {
 } ColumnInfoRecord;
 
 static ColumnInfoRecord sAtomColumns[] = {
-{"atom", 4, 0}, {"name", 4, 1}, {"type", 4, 1}, {"element", 4, 1}, {"residue", 6, 1},
-{"x", 6, 1}, {"y", 6, 1}, {"z", 6, 1}, {"charge", 6, 1}, {NULL}
+       {"atom", 4, 0}, {"name", 4, 1}, {"type", 4, 1}, {"element", 4, 1}, {"residue", 6, 1},
+       {"x", 6, 1}, {"y", 6, 1}, {"z", 6, 1}, {"charge", 6, 1}, {NULL}
 };
 static ColumnInfoRecord sBondColumns[] = {
-{"atoms", 9, 0}, {"names", 9, 0}, {"type", 9, 0}, {"length", 8, 0}, {"r0", 8, 0}, {"force", 8, 0}, {NULL}
+       {"atoms", 9, 0}, {"names", 9, 0}, {"type", 9, 0}, {"length", 8, 0}, 
+       {"par type", 8, 0}, {"k", 8, 0}, {"r0", 8, 0}, {NULL}
 };
 static ColumnInfoRecord sAngleColumns[] = {
-{"atoms", 12, 0}, {"names", 12, 0}, {"type", 12, 0}, {"angle", 8, 0}, {"a0", 8, 0}, {"force", 8, 0}, {NULL}
+       {"atoms", 12, 0}, {"names", 12, 0}, {"type", 12, 0}, {"angle", 8, 0}, 
+       {"par type", 8, 0}, {"k", 8, 0}, {"a0", 8, 0}, {NULL}
 };
 static ColumnInfoRecord sDihedralColumns[] = {
-{"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"dihedral", 8, 0}, {"force", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
+       {"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"dihedral", 8, 0}, 
+       {"par type", 8, 0}, {"k", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
 };
 static ColumnInfoRecord sImproperColumns[] = {
-{"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"improper", 8, 0}, {"force", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
+       {"atoms", 15, 0}, {"names", 15, 0}, {"type", 15, 0}, {"improper", 8, 0}, 
+       {"par type", 8, 0}, {"k", 8, 0}, {"period", 4, 0}, {"phi0", 8, 0}, {NULL}
 };
 static ColumnInfoRecord sParameterColumns[] = {
-{"class", 5, 0}, {"type", 9, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"src", 8, 0}, {"comment", 25, 0}, {NULL}
+       {"class", 5, 0}, {"type", 9, 0}, {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, 
+       {"", 6, 0}, {"", 6, 0}, {"", 6, 0}, {"src", 8, 0}, {"comment", 25, 0}, {NULL}
 };
-static ColumnInfoRecord sMOInfoColumns[] = {
-{"MO", 10, 0}, {"energy", 20, 0}, {NULL}
+static ColumnInfoRecord sUnitCellColumns[] = {
+       {"name", 6, 0}, {"values", 9, 1}, {"", 9, 1}, {"", 9, 1}, {NULL}
+};
+static ColumnInfoRecord sXtalCoordColumns[] = {
+       {"atom", 4, 0}, {"name", 4, 1}, {"element", 4, 1},
+       {"frx", 6, 1}, {"fry", 6, 1}, {"frz", 6, 1},
+       {"occ", 6, 1}, {"temp", 6, 1}, {NULL}
+};
+static ColumnInfoRecord sMOEnergyColumns[] = {
+       {"MO", 5, 0}, {"alpha energy", 12, 0}, {"beta energy", 12, 0}, {NULL}
 };
 static ColumnInfoRecord *sColumnInfo[] = {
-sAtomColumns, sBondColumns, sAngleColumns, sDihedralColumns, sImproperColumns, sParameterColumns, sMOInfoColumns
+       sAtomColumns, sBondColumns, sAngleColumns, sDihedralColumns,
+       sImproperColumns, NULL, sParameterColumns, NULL, 
+       sUnitCellColumns, sXtalCoordColumns, NULL, sMOEnergyColumns
 };
 static char *sTableTitles[] = {
-       "atom", "bond", "angle", "dihedral", "improper", "parameter", "MO info"
+       "atoms", "bonds", "angles", "dihedrals", "impropers", "--------",
+       "MM/MD pars", "--------", "unit cell", "xtal coords", "--------", "MO energy"
 };
 
 void
@@ -2924,7 +3134,7 @@ MainView_tableTitleForIndex(MainView *mview, int idx, char *buf, int bufsize)
        snprintf(buf, bufsize, "%s", sTableTitles[idx]);
 }
 
-void
+int
 MainView_createColumnsForTableAtIndex(MainView *mview, int idx)
 {
        int i;
@@ -2932,8 +3142,10 @@ MainView_createColumnsForTableAtIndex(MainView *mview, int idx)
        if (mview == NULL)
                idx = kMainViewParameterTableIndex;
        if (idx < 0 || idx >= sizeof(sColumnInfo) / sizeof(sColumnInfo[0]))
-               return;
-
+               return 0;  /*  Invalid index  */
+       if (sColumnInfo[idx] == NULL)
+               return 0;  /*  Invalid index  */
+       
        /*  Remove all existing columns  */
        while (MainViewCallback_removeTableColumnAtIndex(mview, 0) > 0);
        
@@ -2951,14 +3163,33 @@ MainView_createColumnsForTableAtIndex(MainView *mview, int idx)
                        width = 80;
                MainViewCallback_addTableColumn(mview, recp->name, width, recp->editable);
        }
+       
+       return 1;
 }
 
 void
 MainView_refreshTable(MainView *mview)
 {
        /*  Reload data  */
-       if (mview != NULL && mview->mol != NULL)
+       if (mview != NULL && mview->mol != NULL) {
                MainView_refreshCachedInfo(mview);
+               if (mview->tableIndex == kMainViewBondTableIndex || 
+                       mview->tableIndex == kMainViewAngleTableIndex || 
+                       mview->tableIndex == kMainViewDihedralTableIndex || 
+                       mview->tableIndex == kMainViewImproperTableIndex || 
+                       mview->tableIndex == kMainViewParameterTableIndex) {
+                       /*  Check the parameter table  */
+                       if (mview->mol->arena == NULL)
+                               md_arena_new(mview->mol);
+                       if (mview->mol->needsMDRebuild || mview->mol->arena->par == NULL) {
+                               /*  Here we do not call MoleculePrepareMDArena(), because we do not want
+                                       to modify Molecule by just redrawing tables.
+                                   (But MoleculePrepareMDArena() *is* called when the parameter
+                                       table is selected.  */
+                               md_prepare(mview->mol->arena, 1);
+                       }
+               }
+       }
        
        MainViewCallback_reloadTableData(mview);
 
@@ -2973,11 +3204,61 @@ MainView_numberOfRowsInTable(MainView *mview)
                return ParameterTableNumberOfRows(gBuiltinParameters);
        if (mview->mol == NULL)
                return 0;
-       if (mview->tableIndex == kMainViewParameterTableIndex)
-               return ParameterTableNumberOfRows(mview->mol->par);
-       if (mview->tableCache == NULL)
-               MainView_refreshCachedInfo(mview);
-       return IntGroupGetCount(mview->tableCache);
+       switch (mview->tableIndex) {
+               case kMainViewParameterTableIndex:
+                       return ParameterTableNumberOfRows(mview->mol->par);
+               case kMainViewUnitCellTableIndex:
+                       return 13; /* a, b, c, alpha, beta, gamma, a_valid, b_valid, c_valid, av, bv, cv, ov */
+               case kMainViewAtomTableIndex:
+               case kMainViewBondTableIndex:
+               case kMainViewAngleTableIndex:
+               case kMainViewDihedralTableIndex:
+               case kMainViewImproperTableIndex:
+               case kMainViewXtalCoordTableIndex:
+                       if (mview->tableCache == NULL)
+                               MainView_refreshCachedInfo(mview);
+                       return IntGroupGetCount(mview->tableCache);
+               case kMainViewMOTableIndex:
+                       return (mview->mol->bset != NULL ? mview->mol->bset->ncomps : 0);
+               default:
+                       return 0;
+       }
+}
+
+int
+MainView_indexToTableRow(MainView *mview, int idx)
+{
+       if (mview == NULL)
+               return -1;
+       switch (mview->tableIndex) {
+               case kMainViewAtomTableIndex:
+               case kMainViewBondTableIndex:
+               case kMainViewAngleTableIndex:
+               case kMainViewDihedralTableIndex:
+               case kMainViewImproperTableIndex:
+               case kMainViewXtalCoordTableIndex:
+                       return IntGroupLookupPoint(mview->tableCache, idx);
+               default:
+                       return idx;
+       }
+}
+
+int
+MainView_tableRowToIndex(MainView *mview, int row)
+{
+       if (mview == NULL)
+               return -1;
+       switch (mview->tableIndex) {
+               case kMainViewAtomTableIndex:
+               case kMainViewBondTableIndex:
+               case kMainViewAngleTableIndex:
+               case kMainViewDihedralTableIndex:
+               case kMainViewImproperTableIndex:
+               case kMainViewXtalCoordTableIndex:
+                       return IntGroupGetNthPoint(mview->tableCache, row);
+               default:
+                       return row;
+       }
 }
 
 static char *
@@ -3044,7 +3325,7 @@ MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufs
                return ParameterTableGetItemText(mview->mol->par, column, row, buf, bufsize);
 
        if (mview->tableIndex == kMainViewAtomTableIndex) { /* Atoms */
-               idx = IntGroupGetNthPoint(mview->tableCache, row);
+               idx = MainView_tableRowToIndex(mview, row);
                ap[0] = ATOM_AT_INDEX(mol->atoms, idx);
                switch (column) {
                        case 0: snprintf(buf, bufsize, "%d", idx); break;
@@ -3064,7 +3345,7 @@ MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufs
                        default: buf[0] = 0; break;
                }
        } else if (mview->tableIndex == kMainViewBondTableIndex) { /* Bonds */
-               idx = IntGroupGetNthPoint(mview->tableCache, row);
+               idx = MainView_tableRowToIndex(mview, row);
                ip = mol->bonds + idx * 2;
                ap[0] = ATOM_AT_INDEX(mol->atoms, ip[0]);
                ap[1] = ATOM_AT_INDEX(mol->atoms, ip[1]);
@@ -3076,17 +3357,23 @@ MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufs
                                snprintf(buf, bufsize, "%.3f", MoleculeMeasureBond(mview->mol, &(ap[0]->r), &(ap[1]->r)));
                                break;
                        case 4:
-                       case 5: {
+                       case 5:
+                       case 6: {
                                BondPar *bp = (BondPar *)sParameterOfTypeAtIndex(mol, kBondParType, idx);
                                if (bp == NULL)
                                        buf[0] = 0;
-                               else snprintf(buf, bufsize, "%.3f", (column == 4 ? bp->r0 : bp->k * INTERNAL2KCAL));
+                               else {
+                                       if (column == 4)
+                                               snprintf(buf, bufsize, "%.6s-%.6s", AtomTypeDecodeToString(bp->type1, typebuf[0]), AtomTypeDecodeToString(bp->type2, typebuf[1]));
+                                       else
+                                               snprintf(buf, bufsize, "%.3f", (column == 5 ? bp->k * INTERNAL2KCAL : bp->r0));
+                               }
                                break;
                        }
                        default: buf[0] = 0; break;
                }
        } else if (mview->tableIndex == kMainViewAngleTableIndex) { /* Angles */
-               idx = IntGroupGetNthPoint(mview->tableCache, row);
+               idx = MainView_tableRowToIndex(mview, row);
                ip = mol->angles + idx * 3;
                ap[0] = ATOM_AT_INDEX(mol->atoms, ip[0]);
                ap[1] = ATOM_AT_INDEX(mol->atoms, ip[1]);
@@ -3099,18 +3386,24 @@ MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufs
                                snprintf(buf, bufsize, "%.3f", MoleculeMeasureAngle(mview->mol, &(ap[0]->r), &(ap[1]->r), &(ap[2]->r)));
                                break;
                        case 4:
-                       case 5: {
+                       case 5:
+                       case 6: {
                                AnglePar *anp = (AnglePar *)sParameterOfTypeAtIndex(mol, kAngleParType, idx);
                                if (anp == NULL)
                                        buf[0] = 0;
-                               else snprintf(buf, bufsize, "%.3f", (column == 4 ? anp->a0 * kRad2Deg : anp->k * INTERNAL2KCAL));
+                               else {
+                                       if (column == 4)
+                                               snprintf(buf, bufsize, "%.6s-%.6s-%.6s", AtomTypeDecodeToString(anp->type1, typebuf[0]), AtomTypeDecodeToString(anp->type2, typebuf[1]), AtomTypeDecodeToString(anp->type3, typebuf[2]));
+                                       else
+                                               snprintf(buf, bufsize, "%.3f", (column == 5 ? anp->k * INTERNAL2KCAL : anp->a0 * kRad2Deg));
+                               }
                                break;
                        }
                        default: buf[0] = 0; break;
                }               
        } else if (mview->tableIndex == kMainViewDihedralTableIndex || mview->tableIndex == kMainViewImproperTableIndex) { /* Dihedrals, Impropers */
                int f = (mview->tableIndex == kMainViewDihedralTableIndex);
-               idx = IntGroupGetNthPoint(mview->tableCache, row);
+               idx = MainView_tableRowToIndex(mview, row);
                ip = (f ? mview->mol->dihedrals : mview->mol->impropers) + idx * 4;
                ap[0] = ATOM_AT_INDEX(mview->mol->atoms, ip[0]);
                ap[1] = ATOM_AT_INDEX(mview->mol->atoms, ip[1]);
@@ -3125,44 +3418,149 @@ MainView_valueForTable(MainView *mview, int column, int row, char *buf, int bufs
                                break;
                        case 4:
                        case 5:
-                       case 6: {
+                       case 6:
+                       case 7: {
                                TorsionPar *tp = (TorsionPar *)sParameterOfTypeAtIndex(mol, (f ? kDihedralParType : kImproperParType), idx);
                                if (tp == NULL)
                                        buf[0] = 0;
-                               else if (column == 5)
+                               else if (column == 4)
+                                       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]));
+                               else if (column == 6)
                                        snprintf(buf, bufsize, "%d", tp->period[0]);
-                               else snprintf(buf, bufsize, "%.3f", (column == 4 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
+                               else snprintf(buf, bufsize, "%.3f", (column == 5 ? tp->k[0] * INTERNAL2KCAL : tp->phi0[0] * kRad2Deg));
                                break;
                        }
                        default: buf[0] = 0; break;
                }               
-       } else if (mview->tableIndex == kMainViewMOTableIndex) { /* MO info */
+       } else if (mview->tableIndex == kMainViewUnitCellTableIndex) { /* Unit cell info */
+               buf[0] = 0;
+               if (mview->mol->cell != NULL) {
+                       static const char *unitCellRowTitles[] = {"a", "b", "c", "alpha", "beta", "gamma", "a_valid", "b_valid", "c_valid", "a_vec", "b_vec", "c_vec", "origin"};
+                       if (column == 0)
+                               snprintf(buf, bufsize, "%s", unitCellRowTitles[row]);
+                       else {
+                               switch(row) {
+                                       case 0: case 1: case 2:
+                                               if (column == 1)
+                                                       snprintf(buf, bufsize, "%.5f", mview->mol->cell->cell[row]);
+                                               else if (column == 2 && mview->mol->cell->has_sigma)
+                                                       snprintf(buf, bufsize, "%.5f", mview->mol->cell->cellsigma[row]);                                                       
+                                               break;
+                                       case 3: case 4: case 5:
+                                               if (column == 1)
+                                                       snprintf(buf, bufsize, "%.4f", mview->mol->cell->cell[row]);
+                                               else if (column == 2 && mview->mol->cell->has_sigma)
+                                                       snprintf(buf, bufsize, "%.4f", mview->mol->cell->cellsigma[row]);                                                       
+                                               break;
+                                       case 6: case 7: case 8:
+                                               if (column == 1)
+                                                       snprintf(buf, bufsize, "%d", (int)mview->mol->cell->flags[row - 6]);
+                                               break;
+                                       case 9: case 10: case 11: case 12: {
+                                               Vector *vp = (row == 12 ? &mview->mol->cell->origin : &mview->mol->cell->axes[row - 9]);
+                                               Double dval = (column == 1 ? vp->x : (column == 2 ? vp->y : vp->z));
+                                               snprintf(buf, bufsize, "%.5f", dval);
+                                               break;
+                                       }
+                               }
+                       }
+               }
+       } else if (mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* fractional coords */
+               Vector fract_r;
+               idx = MainView_tableRowToIndex(mview, row);
+               ap[0] = ATOM_AT_INDEX(mol->atoms, idx);
+               if (mview->mol->cell != NULL)
+                       TransformVec(&fract_r, mview->mol->cell->rtr, &(ap[0]->r));
+               else fract_r = ap[0]->r;
+               switch (column) {
+                       case 0: snprintf(buf, bufsize, "%d", idx); break;
+                       case 1: snprintf(buf, bufsize, "%.4s", ap[0]->aname); break;
+                       case 2: snprintf(buf, bufsize, "%.2s", ap[0]->element); break;
+                       case 3: snprintf(buf, bufsize, "%.3f", fract_r.x); break;
+                       case 4: snprintf(buf, bufsize, "%.3f", fract_r.y); break;
+                       case 5: snprintf(buf, bufsize, "%.3f", fract_r.z); break;
+                       case 6: snprintf(buf, bufsize, "%.3f", ap[0]->occupancy); break;
+                       case 7: snprintf(buf, bufsize, "%.3f", ap[0]->tempFactor); break;
+                       default: buf[0] = 0; break;
+               }               
+       } else if (mview->tableIndex == kMainViewMOTableIndex) { /* MO energy */
                BasisSet *bset = mview->mol->bset;
                idx = row;
                buf[0] = 0;
-               if (bset != NULL && idx >= 0 && idx < bset->nmos) {
-                       char *s = "";
-                       if (idx * 2 + 2 == bset->nelectrons)
-                               s = " (HOMO)";
-                       else if (idx * 2 + 1 == bset->nelectrons)
-                               s = " (SOMO)";
-                       else if ((bset->nelectrons + 1) / 2 == idx)
-                               s = " (LUMO)";
-                       switch (column) {
-                               case 0: snprintf(buf, bufsize, "%d%s", idx + 1, s); break;
-                               case 1: snprintf(buf, bufsize, "%.8f", bset->moenergies[idx]); break;
+               if (bset != NULL && idx >= 0 && idx < bset->ncomps) {
+                       if (column == 0) {
+                               snprintf(buf, bufsize, "%d", idx + 1);
+                       } else if (column >= 3) {
+                               return;  /*  Cannot happen  */
+                       } else {
+                               /*  column == 1, alpha; column == 2, beta  */
+                               char eflag = ' ';
+                               if (column == 1) {
+                                       /*  Alpha  */
+                                       if (idx >= bset->ne_alpha)
+                                               eflag = '*';  /*  Unoccupied  */
+                                       else {
+                                               if (bset->rflag == 2 && idx >= bset->ne_beta)
+                                                       eflag = 'S';  /*  Singly occupied  */
+                                       }
+                               } else {
+                                       /*  Beta  */
+                                       if (bset->rflag != 0)
+                                               return;  /*  No beta orbitals  */
+                                       if (idx >= bset->ne_beta)
+                                               eflag = '*';  /*  Unoccupied  */
+                                       idx += bset->ncomps;
+                               }
+                               snprintf(buf, bufsize, "%c %.8f", eflag, bset->moenergies[idx]);
                        }
                }
        }
 }
 
-/*  Set color for the locally defined or undefined MM parameters  */
+/*  Set color for the table  */
 int
 MainView_setColorForTable(MainView *mview, int column, int row, float *fg, float *bg)
 {
        int parType = -1;
        int idx;
        UnionPar *up;
+
+       if (mview->tableIndex == kMainViewParameterTableIndex && column == -1) {
+               int src = ParameterTableGetItemSource(mview->mol->par, row);
+               if (src == -2) {  /* separator line */
+                       bg[0] = bg[1] = bg[2] = 0.6;
+                       return 2;
+               } else if (src == -1) { /*  undefined parameters  */
+                       bg[0] = 1.0;
+                       bg[1] = bg[2] = 0.2;
+                       return 2;
+               } else if (src == 0) {  /*  local parameter  */
+                       bg[0] = bg[1] = 1.0;
+                       bg[2] = 0.6;
+                       return 2;
+               }
+       } else if (mview->tableIndex == kMainViewMOTableIndex && column == -1) {
+               BasisSet *bset = mview->mol->bset;
+               int n = 0;
+               if (bset == NULL)
+                       return 0;
+               if (row < 0 || row >= bset->ncomps)
+                       return 0;
+               if (row < bset->ne_beta)
+                       n = 0;  /*  Occupied  */
+               else if (row < bset->ne_alpha)
+                       n = 1;  /*  Half-occupied  */
+               else n = 2;  /*  Unoccupied  */
+               switch (n) {
+                       case 1: bg[0] = bg[1] = bg[2] = 0.9; break;
+                       case 2: bg[0] = bg[1] = bg[2] = 0.8; break;
+                       default: bg[0] = bg[1] = bg[2] = 1.0; break;
+               }
+               return 2;
+       } else if (mview->tableIndex < kMainViewBondTableIndex || mview->tableIndex > kMainViewImproperTableIndex)
+               return 0;
+       
+       /*  Bond etc. table  */
        switch (mview->tableIndex) {
                case kMainViewBondTableIndex:
                        parType = kBondParType;
@@ -3182,7 +3580,7 @@ MainView_setColorForTable(MainView *mview, int column, int row, float *fg, float
        if (column < 3)
                return 0;
 
-       idx = IntGroupGetNthPoint(mview->tableCache, row);
+       idx = MainView_tableRowToIndex(mview, row);
        up = sParameterOfTypeAtIndex(mview->mol, parType, idx);
        if (up == NULL)
                return 0;
@@ -3222,8 +3620,17 @@ MainView_setSelectionFromTable(MainView *mview)
        IntGroup *ig, *sel;
        if (mview == NULL || mview->mol == NULL)
                return;
-       if (mview->tableIndex == kMainViewMOTableIndex || mview->tableIndex == kMainViewParameterTableIndex)
-               return;  /*  Do nothing  */
+       switch (mview->tableIndex) {
+               case kMainViewAtomTableIndex:
+               case kMainViewBondTableIndex:
+               case kMainViewAngleTableIndex:
+               case kMainViewDihedralTableIndex:
+               case kMainViewImproperTableIndex:
+               case kMainViewXtalCoordTableIndex:
+                       break;   /*  Continue  */
+               default:
+                       return;  /*  Do nothing  */
+       }
        ig = MainViewCallback_getTableSelection(mview);
        sel = IntGroupNew();
        if (ig != NULL) {
@@ -3233,7 +3640,8 @@ MainView_setSelectionFromTable(MainView *mview)
                        i2 = IntGroupGetNthPoint(mview->tableCache, i1);
                        if (i2 < 0)
                                continue;
-                       if (mview->tableIndex == kMainViewAtomTableIndex) {  /* Atoms */
+                       if (mview->tableIndex == kMainViewAtomTableIndex ||
+                               mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Atoms */
                                IntGroupAdd(sel, i2, 1);
                        } else if (mview->tableIndex == kMainViewBondTableIndex) {  /* Bonds */
                                ip = mview->mol->bonds + i2 * 2;
@@ -3335,7 +3743,7 @@ MainView_setValueForTable(MainView *mview, int column, int row, const char *buf)
        MainView_valueForTable(mview, column, row, temp, sizeof temp);
        if (strcmp(buf, temp) == 0 || buf[0] == 0)
                return;  /*  No change  */
-       idx = IntGroupGetNthPoint(mview->tableCache, row);
+       idx = MainView_tableRowToIndex(mview, row);
        if (mview->tableIndex == kMainViewAtomTableIndex) { /* Atoms */
                switch (column) {
                        case 1: key = "name"; break;
@@ -3356,6 +3764,62 @@ MainView_setValueForTable(MainView *mview, int column, int row, const char *buf)
                MolActionCreateAndPerform(mol, SCRIPT_ACTION("iss"), "set_atom_attr", idx, key, buf);
        } else if (mview->tableIndex == kMainViewParameterTableIndex) { /* Parameters */
                sMainView_ParameterTableSetItemText(mview->mol, column, row, buf);
+       } else if (mview->tableIndex == kMainViewUnitCellTableIndex) { /* Unit cell */
+               Double cellPars[12], dval;
+               char flags[3];
+               int n1;
+               Vector vecs[4];
+               if (mol->cell == NULL) {
+                       /*  Create a cell  */
+                       static const Double sCellPars[] = {1, 1, 1, 90, 90, 90};
+                       MolActionCreateAndPerform(mol, gMolActionSetCell, 6, sCellPars, 0);
+               }
+               if (row >= 0 && row < 6) {
+                       n1 = 6;
+                       memmove(cellPars, mol->cell->cell, sizeof(Double) * 6);
+                       if (mol->cell->has_sigma)
+                               memmove(cellPars + 6, mol->cell->cellsigma, sizeof(Double) * 6);
+                       else memset(cellPars + 6, 0, sizeof(Double) * 6);
+                       dval = strtod(buf, NULL);
+                       if (column == 1) {
+                               if (dval == 0.0)
+                                       return;
+                               cellPars[row] = dval;
+                       } else {
+                               n1 = 12;
+                               cellPars[row + 6] = dval;
+                       }
+                       MolActionCreateAndPerform(mol, gMolActionSetCell, n1, cellPars, 0);
+               } else {
+                       memmove(vecs, mol->cell->axes, sizeof(Vector) * 3);
+                       vecs[3] = mol->cell->origin;
+                       memmove(flags, mol->cell->flags, 3);
+                       if (row >= 6 && row < 9)
+                               flags[row - 6] = strtol(buf, NULL, 0);
+                       else {
+                               Vector *vp = vecs + (row - 9);
+                               dval = strtod(buf, NULL);
+                               if (column == 1)
+                                       vp->x = dval;
+                               else if (column == 2)
+                                       vp->y = dval;
+                               else vp->z = dval;
+                       }
+                       MolActionCreateAndPerform(mol, gMolActionSetBox, vecs, vecs + 1, vecs + 2, vecs + 3, (flags[0] != 0) * 4 + (flags[1] != 0) * 2 + (flags[2] != 0), 0);
+               }
+               
+       } else if (mview->tableIndex == kMainViewXtalCoordTableIndex) {  /* Fractional coords */
+               switch (column) {
+                       case 1: key = "name"; break;
+                       case 2: key = "element"; break;
+                       case 3: key = "fract_x"; break;
+                       case 4: key = "fract_y"; break;
+                       case 5: key = "fract_z"; break;
+                       case 6: key = "occupancy"; break;
+                       case 7: key = "temp_factor"; break;
+                       default: return;
+               }
+               MolActionCreateAndPerform(mol, SCRIPT_ACTION("iss"), "set_atom_attr", idx, key, buf);
        }
 }
 
@@ -3368,6 +3832,13 @@ MainView_isTableItemEditable(MainView *mview, int column, int row)
                return 0;
        if (mview->tableIndex == kMainViewParameterTableIndex)
                return ParameterTableIsItemEditable(mview->mol->par, column, row);
+       if (mview->tableIndex == kMainViewUnitCellTableIndex) {
+               if ((row >= 0 && row < 6 && column >= 1 && column <= 3) ||
+                       (row >= 6 && row < 9 && column == 1) ||
+                       (row >= 9 && row < 13 && column >= 1 && column <= 3))
+                       return 1;
+               else return 0;
+       }
        if (mview->tableIndex >= 0 && mview->tableIndex < sizeof(sColumnInfo) / sizeof(sColumnInfo[0]))
                return sColumnInfo[mview->tableIndex][column].editable != 0;
        else return 0;
@@ -3386,14 +3857,24 @@ MainView_tableType(MainView *mview)
 void
 MainView_dragTableSelectionToRow(MainView *mview, int row)
 {
-       Int *new2old, i, n, count, natoms, start_row;
-       IntGroup *sel;
-       if (mview == NULL || mview->mol == NULL || mview->tableIndex != 0 || row < 0 || row > (natoms = mview->mol->natoms))
+       Int *new2old, i, n, count, natoms, start_idx, to_idx;
+    IntGroup *sel, *sel2;
+       if (mview == NULL || mview->mol == NULL)
                return;
+       if (mview->tableIndex != kMainViewAtomTableIndex && mview->tableIndex != kMainViewXtalCoordTableIndex)
+               return;  /*  Only atom tables can respond  */
        if (md_is_running(mview->mol->arena)) {
                MoleculeCallback_cannotModifyMoleculeDuringMDError(mview->mol);
                return;
-       }       
+       }
+       n = MainView_numberOfRowsInTable(mview);
+       if (row < 0 || row > n)
+               return;  /*  Out of range  */
+       natoms = mview->mol->natoms;
+       if (row == n)
+               to_idx = natoms;
+       else
+               to_idx = MainView_tableRowToIndex(mview, row);
        sel = MoleculeGetSelection(mview->mol);
        if (sel == NULL || (count = IntGroupGetCount(sel)) == 0)
                return;
@@ -3402,34 +3883,56 @@ MainView_dragTableSelectionToRow(MainView *mview, int row)
                return;
 
        //  Place the atoms above the target position
-       for (i = n = 0; i < row; i++) {
+       for (i = n = 0; i < to_idx; i++) {
                if (IntGroupLookupPoint(sel, i) < 0)
                        new2old[n++] = i;
        }
-       start_row = n;
+       start_idx = n;
        //  Place the atoms within the selection
        for (i = 0; i < count; i++) {
                new2old[n++] = IntGroupGetNthPoint(sel, i);
        }
        //  Place the remaining atoms
-       for (i = row; i < natoms; i++) {
+       for (i = to_idx; i < natoms; i++) {
                if (IntGroupLookupPoint(sel, i) < 0)
                        new2old[n++] = i;
        }
        MolActionCreateAndPerform(mview->mol, gMolActionRenumberAtoms, n, new2old);
        
        //  Change selection
-       sel = IntGroupNewWithPoints(start_row, count, -1);
+    //  Molecule selection (atom indices)
+       sel = IntGroupNewWithPoints(start_idx, count, -1);
        MolActionCreateAndPerform(mview->mol, gMolActionSetSelection, sel);
+    //  Table selection (row numbers)
+    sel2 = IntGroupNew();
+    for (i = 0; i < count; i++) {
+        int row_i = MainView_indexToTableRow(mview, i + start_idx);
+        if (row_i >= 0)
+            IntGroupAdd(sel2, row_i, 1);
+    }
+    MainViewCallback_setTableSelection(mview, sel2);
+    IntGroupRelease(sel2);
        IntGroupRelease(sel);
 }
 
+int
+MainView_isRowSelectable(MainView *mview, int row)
+{
+    if (mview->tableIndex == kMainViewParameterTableIndex) {
+        int src = ParameterTableGetItemSource(mview->mol->par, row);
+        if (src == -2) {  /* separator line */
+            return 0;
+        }
+    }
+    return 1;
+}
+
 IntGroup *
 MainView_selectedMO(MainView *mview)
 {
+    if (!gUseGUI)
+        return NULL;
        if (mview == NULL || mview->mol == NULL || mview->tableIndex != kMainViewMOTableIndex)
                return NULL;
        return MainViewCallback_getTableSelection(mview);  /*  Note: the indices are 0 based  */
 }
-
-       
\ No newline at end of file