/*
 ----------------------------------------------------------------------------
 "THE BEER-WARE LICENSE" (Revision 42):
 <daniel.kratzert@ac.uni-freiburg.de> wrote this file. As long as you retain
 this notice you can do whatever you want with this stuff. If we meet some day,
 and you think this stuff is worth it, you can buy me a beer in return.
 Daniel Kratzert
 ----------------------------------------------------------------------------
*/

#include "dsrglwindow.h"
#include "molecule.h"
#include "chgl.h"
#include "dsrgui.h"


DSRGlWindow::DSRGlWindow(QWidget *parent, Molecule *m,
                         DSRMol header, QString fragment) {
  this->setParent(parent);
  m_molecule=m;
  showFit=NULL;
  showFitLabel=NULL;
  mole.cell.symmops.clear();
  mole.cell.trans.clear();
  mole.asymm.clear();
  mole.sdm.clear();
  mole.showbonds.clear();
  mole.showatoms.clear();
  mole.selectedatoms.clear();
  mole.envi_sdm.clear();
  mole.contact.clear();
  mole.envi_sdm.clear();
  mole.contact.clear();
  mole.symmopsEQIV.clear();
  mole.labelEQIV.clear();
  mole.transEQIV.clear();
  mole.freeatoms.clear();
  mole.bindatoms.clear();
  mole.selectedatoms.clear();
  mole.showatoms.clear();
  dsrLabelColor = QColor("white");
  V3 nl(0,0,0);
  mole.cell.trans.append(nl);
  mole.cell.symmops.append(Matrix(1,0,0, 0,1,0, 0,0,1));
  mole.pmin=10000;
  mole.pmax=-10000;
  mole.LOD=3;
  // initial zoom level:
  // The Zom will later be adapted to fragment size.
  L=50.0;
  mole.transEQIV.clear();
  mole.symmopsEQIV.clear();
  mole.labelEQIV.clear();
  mole.freeatoms.clear();
  mole.bindatoms.clear();
  mole.envi_sdm.clear();
  if (fragment.size() > 0) {
    display_fragment(header);
  }
}

DSRGlWindow::~DSRGlWindow() {
}

static inline void transform_point(GLdouble out[4], const GLdouble m[16], const GLdouble in[4]) {
#define M(row,col)  m[col*4+row]
  out[0] =
	  M(0, 0) * in[0] + M(0, 1) * in[1] + M(0, 2) * in[2] + M(0, 3) * in[3];
  out[1] =
	  M(1, 0) * in[0] + M(1, 1) * in[1] + M(1, 2) * in[2] + M(1, 3) * in[3];
  out[2] =
	  M(2, 0) * in[0] + M(2, 1) * in[1] + M(2, 2) * in[2] + M(2, 3) * in[3];
  out[3] =
	  M(3, 0) * in[0] + M(3, 1) * in[1] + M(3, 2) * in[2] + M(3, 3) * in[3];
#undef M
}

static inline bool  posTo2D(V3 obj,
		const GLdouble model[16], const GLdouble proj[16],
		const GLint viewport[4],
		GLdouble * winx, GLdouble * winy) {
  GLdouble in[4], out[4];

  in[0] = obj.x;
  in[1] = obj.y;
  in[2] = obj.z;
  in[3] = 1.0;
  transform_point(out, model, in);
  transform_point(in, proj, out);

  if (in[3] == 0.0) return false;

  in[0] /= in[3];
  in[1] /= in[3];
  in[2] /= in[3];

  *winx = viewport[0] + (1 + in[0]) * viewport[2] / 2;
  *winy = viewport[1] + (1 - in[1]) * viewport[3] / 2;
  return true;
}


void DSRGlWindow::draw(){
  /*const GLfloat  OBJ_SPE[]   = { 0.8, 0.8, 0.8, 1.0 };
    const GLfloat  OBJ_SHIN    = 127.0;               
    glMaterialfv( GL_FRONT_AND_BACK, GL_SPECULAR,             OBJ_SPE  );
    glEnable     ( GL_COLOR_MATERIAL ) ;
    glColorMaterial ( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE ) ;    
    glMaterialf(  GL_FRONT_AND_BACK, GL_SHININESS,           OBJ_SHIN );
  */
  mole.dratom = 0;
  glPushMatrix();
  glScaled( L, L, L );
  GLdouble model[16];
  GLdouble proj[16];
  GLint viewport[4];
  glGetIntegerv(GL_VIEWPORT, viewport);
#if defined Q_OS_MAC && (QT_VERSION >= 0x050000)
  viewport[2]/=2;
  viewport[3]/=2;
#endif
  //      printf("%d %d %d %d\n",viewport[0],viewport[1],viewport[2],viewport[3]);
  glGetDoublev( GL_PROJECTION_MATRIX, (double*)proj );
  glGetDoublev( GL_MODELVIEW_MATRIX, (double*)model );
  mole.adp = 0;
  mole.intern = 0;
  mole.tubes = 0;
  glDisable(GL_BLEND);
  mole.atoms(xd, 50);
  mole.bonds(bonds);
  glEnable(GL_COLOR_MATERIAL);
  qglColor(dsrLabelColor);  // The fragment label color
  glClear( GL_DEPTH_BUFFER_BIT);
  for (int i=0; i<xd.size(); i++){
      GLdouble in[4], out[4];
      in[0] = xd.at(i).pos.x;
      in[1] = xd.at(i).pos.y;
      in[2] = xd.at(i).pos.z;
      in[3] = 1.0;
      transform_point(out, model, in);
      posTo2D(xd.at(i).pos,model,proj,viewport, &xd[i].screenX, &xd[i].screenY);
      bool issel=false;
      for (int j=0; j<selFragAt.size(); j++) {
        if ((showFitLabel!=NULL)
              && (!showFitLabel->isChecked())
              && (j<m_molecule->selectedatoms.size())
              && (xd.at(i).Label==selFragAt.at(j).Label)) {
          issel=true;
          renderText( xd.at(i).pos.x, xd.at(i).pos.y, xd.at(i).pos.z,
              xd.at(i).Label+"->"+m_molecule->selectedatoms.at(j).Label, myFont);
        }
      }
      if (!issel) {
        renderText( xd.at(i).pos.x, xd.at(i).pos.y, xd.at(i).pos.z,
            xd.at(i).Label, myFont);
      }
  }
  glPopMatrix();
  if (!selFragAt.isEmpty()){
    glPushMatrix();{
      glScaled( L, L, L );
      mole.tubes=0;
      mole.intern=1;
      mole.adp=0;
      mole.dratom=1;
      mole.atoms(selFragAt);
      mole.dratom=0;
      glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);
    }glPopMatrix();
  }
  if ((showFit!=NULL)&&(!fita.isEmpty())&&(showFit->isChecked())){
    glPushMatrix();{
      glEnable(GL_BLEND);
      glScaled( L, L, L );
      mole.tubes=0;
      mole.intern=0;
      mole.adp=0;
      mole.dratom=5;
      mole.atoms(fita);
      mole.dbond(fitabonds);
      mole.lbond();
      if (showFitLabel->isChecked()){
        qglColor(QColor(0, 0, 175));  // The surrounding atoms label color
        for (int i=0; i<fita.size(); i++){
          renderText( fita.at(i).pos.x, fita.at(i).pos.y, fita.at(i).pos.z, fita.at(i).Label, myFont);
        }
      }
    }glPopMatrix();
  }
}

void DSRGlWindow::initializeGL(){
  glEnable(GL_LINE_SMOOTH);  
  glHint(GL_LINE_SMOOTH_HINT,GL_NICEST);
  //glEnable(GL_POLYGON_SMOOTH);   
  myFont = QFont("Arial", 12, -1, false);  // the initial label font size
  const GLfloat  position[] = {100.0f, 100.0f,100.0f,0.0f};
  const GLfloat  diffuse[]  = { 1.0, 1.0, 1.0, 1.0 };
  const GLfloat  specular[] = { 1.0, 0.9, 0.9, 1.0 };
  const GLfloat  ambient[]  = { 0.4, 0.4, 0.4, 1.0 };
  glLightModeli(  GL_LIGHT_MODEL_LOCAL_VIEWER, 1 );
  glLightfv( GL_LIGHT0, GL_POSITION, position );
  glLightfv( GL_LIGHT0, GL_AMBIENT,  ambient );
  glLightfv( GL_LIGHT0, GL_DIFFUSE,  diffuse );
  glLightfv( GL_LIGHT0, GL_SPECULAR, specular );
  glLightfv( GL_LIGHT1, GL_POSITION, position );
  glLightfv( GL_LIGHT1, GL_DIFFUSE,  diffuse );  
  glLightfv( GL_LIGHT1, GL_AMBIENT,  ambient );
  glLightfv( GL_LIGHT2, GL_DIFFUSE,  diffuse );  
  glEnable( GL_LIGHTING );
  glEnable( GL_LIGHT0 );
  //  glEnable( GL_BLEND);
  glDisable(GL_BLEND);
  glAlphaFunc ( GL_GREATER, 0.01f ) ;
  //glEnable(GL_ALPHA_TEST);
  glBlendFunc ( GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA ) ;
  const GLfloat  OBJ_SPE[]   = { 1.0, 1.0, 1.0, 1.0 };  
  const GLfloat  OBJ_SHIN    = 127.0;                   
  glMaterialfv( GL_FRONT_AND_BACK, GL_SPECULAR,             OBJ_SPE  );
  glEnable     ( GL_COLOR_MATERIAL ) ;
  glColorMaterial ( GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE ) ;    
  glMaterialf(  GL_FRONT_AND_BACK, GL_SHININESS,           OBJ_SHIN );
  glShadeModel( GL_SMOOTH );
  glEnable(GL_NORMALIZE);
  glClearColor(0.6,0.6,0.6,1.0);
  glEnable(GL_DEPTH_TEST );
  glDepthFunc(GL_LEQUAL);
  gluLookAt_(0.0, 200, 50 ,   0.0, 0.0, 0.0,   0.0, -100, 400 );
  glHint(GL_PERSPECTIVE_CORRECTION_HINT, GL_NICEST);
}

void DSRGlWindow::resizeGL(int width, int height){
#if defined Q_OS_MAC && (QT_VERSION >= 0x050000)
  glViewport(0, 0, width*2,height*2);        
#else  
  glViewport(0, 0, width, height);        
#endif
  glGetIntegerv(GL_VIEWPORT, vp);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective_( 29.0, (double)width/height, 5.0, 8000.0 );
}

QSize DSRGlWindow::minimumSizeHint() const
{
    return QSize(250, 250);
}

QSize DSRGlWindow::sizeHint() const
{
    return QSize(400, 400);
}

void DSRGlWindow::paintGL(){
  glClearColor(0.6,0.6,0.6,1.0); // The background color
  glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
#if defined Q_OS_MAC && (QT_VERSION >= 0x050000)
  glViewport(0, 0, QGLWidget::width()*2, QGLWidget::height()*2);        
#else  
  glViewport(0, 0, QGLWidget::width(), QGLWidget::height());        
#endif
  glGetIntegerv(GL_VIEWPORT, vp);
  glMatrixMode(GL_PROJECTION);
  glLoadIdentity();
  gluPerspective_( 29.0, (double)QGLWidget::width()/QGLWidget::height(), 5.0, 8000.0 );
  glMatrixMode(GL_MODELVIEW);
  glPushMatrix();
  draw();
  glPopMatrix();
}

void DSRGlWindow::clear_molecule()
{ // clear the 3D view
  xd.clear();
  bonds.clear();
  fita.clear();
  fitabonds.clear();
  selFragAt.clear();
  source_atoms.clear();
  updateGL();
}

void DSRGlWindow::display_fragment(DSRMol frag, bool clear_sel)
{ //! uses a QStringlist of atoms from DSR and fills
  //! them into a MyAtom. MyAtom is then drawn in OpenGL
  //! @param The fragment from DSR
  //! @param clear_sel if enabled, the selected atoms are cleared
  //! before drawing the molecule
  V3 mid = V3(0,0,0);
  MyAtom newAtom;
  newAtom.part = 0;
  newAtom.resiNr = 0;
  newAtom.hidden = 0;
  newAtom.symmGroup = 0;
  newAtom.afix = 0;
  xd.clear();
  bonds.clear();
  fita.clear();
  fitabonds.clear();
  // Do not clear the selection during invert:
  if (clear_sel) {
    selFragAt.clear();
  }
  if (frag.atoms.isEmpty()) {
    return;
  }
  foreach (QStringList atom, frag.atoms) {
    if (atom.length() < 4)  continue;  // Do not crash on empty string
    newAtom.Label = atom.at(0);
    newAtom.an = abs(atom.at(1).toInt())-1;
    newAtom.pos.x = atom.at(2).toDouble();
    newAtom.pos.y = atom.at(3).toDouble();
    newAtom.pos.z = atom.at(4).toDouble();
    mid += newAtom.pos;
    //qDebug()<< newAtom.Label << newAtom.an << newAtom.pos.x << newAtom.pos.y << newAtom.pos.z;
    xd.append(newAtom);
  }
  mid *= -1.0/xd.size();
  for (int i=0; i<xd.size(); i++){
    xd[i].pos += mid;
  }
  L = 90.0/mole.dimension(xd); // Zoom level of the molecule inside the widget
  bonds=mole.connecting(xd,true);
  updateGL();
}


void DSRGlWindow::set_label_color() {
    //! Sets the label color. The default color was not visible
    //! for white or yellow atoms.
    QColor lc= QColorDialog::getColor(QColor(6, 167, 222), this);
    if (lc.isValid()) {
        dsrLabelColor=lc;
    } else {
        dsrLabelColor = QColor("white");
    }
    updateGL();
}


void DSRGlWindow::mousePressEvent(QMouseEvent *event){
  //! Handles the mousePressEvents in the 3D view.
  //! Atoms are selectable if showFit is not NULL.
  //! New atoms will be selected if the mousePressEvent is inside the
  //! size of the respective atom and less than three atoms.
  //! With more than three atoms, the first selection will be removed.
  moux=event->pos().x();
  mouy=event->pos().y();
  if((event->buttons() == Qt::RightButton && event->modifiers() == Qt::ControlModifier)){
      QMenu *menu = new QMenu("");
      menu->addAction(QString("Set label color"), this, SLOT(set_label_color()));
      menu->exec(QCursor::pos());
  }
  if (event->buttons() == Qt::LeftButton) {
    // Has to be after moux/mouy, otherwise zoom is strange:
    if (showFit == NULL) return;  // <-prevents selection of atoms in edit window
    double nahda=200.0, da=0;
    int nahdai=-1;
    for (int j=0; j<xd.size();j++){
      if (xd.at(j).hidden) continue;
      da=(((xd.at(j).screenX-event->x())*( xd.at(j).screenX-event->x()))+
          ((xd.at(j).screenY-event->y())*( xd.at(j).screenY-event->y())));
      nahdai=(da<nahda)?j:nahdai;
      nahda=qMin(nahda,da);
    }
    if ((nahdai >- 1) && (nahdai < xd.size())){
      // This might be a good idea in the first place, but results in wired behavior:
      //if (selFragAt.contains(xd[nahdai])) {
      //    return;
      //}
      selFragAt.append(xd[nahdai]);
      if (selFragAt.size() > 3) selFragAt.removeFirst();
      updateGL();
      if (selFragAt.size() == 3) {
        source_atoms = QString("%1 %2 %3")
                        .arg(selFragAt.at(0).Label)
                        .arg(selFragAt.at(1).Label)
                        .arg(selFragAt.at(2).Label);
        emit sourceStringChanged();
      }
    }
    makeInfo();
  }
}

void DSRGlWindow::makeInfo(){
  //! makes an info box that shows the fragment fit distances
  double f12dis=0, f13dis=0, t12dis=0, t13dis=0;
  if (selFragAt.size()>1) f12dis=sqrt(Distance(selFragAt.at(0).pos,selFragAt.at(1).pos));
  if (selFragAt.size()>2) f13dis=sqrt(Distance(selFragAt.at(0).pos,selFragAt.at(2).pos));
  if (m_molecule->selectedatoms.size()>1) t12dis=sqrt(Distance(m_molecule->selectedatoms.at(0).pos,m_molecule->selectedatoms.at(1).pos));
  if (m_molecule->selectedatoms.size()>2) t13dis=sqrt(Distance(m_molecule->selectedatoms.at(0).pos,m_molecule->selectedatoms.at(2).pos));
  emit updateInfo(QString(
        "<table>"
        "<tr><th colspan=\"2\" align=\"left\">fragment: </th>"
        "<th><th>"
        "    <th colspan=\"2\" align=\"left\">target:</th>"
        "</tr>"
        "<tr><td colspan=\"2\"><font color=blue>%1</font></td>"
        "<td><td>"
        "    <td colspan=\"2\"><font color=blue>%2</font></td>"
        "</tr>"
        "<tr><td align=\"left\"><font color=blue>%3</font></td> "
        "    <td align=\"left\">1,2: <b>%7 &Aring;</b></td>"
        "<td><td>"
        "    <td align=\"left\"><font color=blue>%4 </font></td>"
        "    <td align=\"right\">1,2: <b>%9 &Aring;</b></td>"
        "</tr>"
        "<tr><td align=\"left\"><font color=blue>%5</font></td>"
        "    <td align=\"left\">1,3: <b>%8 &Aring;</b></td>"
        "<td><td>"
        "    <td align=\"right\"><font color=blue>%6 </font></td>"
        "    <td align=\"right\">1,3: <b>%10 &Aring;<b></td>"
        "</tr>"
        "</table>")
    .arg((selFragAt.size()>0)?selFragAt.at(0).Label:"")
    .arg((m_molecule->selectedatoms.size()>0)?m_molecule->selectedatoms.at(0).Label:"")
    .arg((selFragAt.size()>1)?selFragAt.at(1).Label:"")
    .arg((m_molecule->selectedatoms.size()>1)?m_molecule->selectedatoms.at(1).Label:"")
    .arg((selFragAt.size()>2)?selFragAt.at(2).Label:"")
    .arg((m_molecule->selectedatoms.size()>2)?m_molecule->selectedatoms.at(2).Label:"")
    .arg(f12dis, 4, 'f', 2)
    .arg(f13dis, 4, 'f', 2)
    .arg(t12dis, 4, 'f', 2)
    .arg(t13dis, 4, 'f', 2)
    );
    tryfit();
}

//void glRotateL(double ang, double x, double y, double z);
inline void DSRGlWindow::__RotateCS( double c, double s, double& X, double& Y ) {
  double T = X;
  X = c*X - s*Y;
  Y = s*T + c*Y;
}

void DSRGlWindow::glTranslateL( const double dx, const double dy, const double dz ) {
  double mat[4][4];

  glGetDoublev( GL_MODELVIEW_MATRIX, (double*)mat );
  mat[3][0] += dx;  mat[3][1] += dy;  mat[3][2] += dz;
  glLoadMatrixd((double*)mat);
}
void DSRGlWindow::glRotateL( const double dang, const double x, const double y, const double z ) {
  double mat[4][4];
#ifndef M_PI
#define	M_PI		3.14159265358979323846	/* pi */
#endif

  double s = z;
  s = sin(dang*M_PI/180);
  const double c = cos(dang*M_PI/180);
  glGetDoublev( GL_MODELVIEW_MATRIX, (double*)mat );
  //  glGetDoublev( GL_PROJECTION_MATRIX, (double*)mat );
  if( x!=0.0 ){
    __RotateCS( c, s, mat[0][1], mat[0][2] );
    __RotateCS( c, s, mat[1][1], mat[1][2] );
    __RotateCS( c, s, mat[2][1], mat[2][2] );
    //    printf("x\n");
  }else if( y!=0.0 ){
    __RotateCS( c, s, mat[0][2], mat[0][0] );
    __RotateCS( c, s, mat[1][2], mat[1][0] );
    __RotateCS( c, s, mat[2][2], mat[2][0] );
    //    printf("y\n");
  }else{
    __RotateCS( c, s, mat[0][0], mat[0][1] );
    __RotateCS( c, s, mat[1][0], mat[1][1] );
    __RotateCS( c, s, mat[2][0], mat[2][1] );
    //    printf("z\n");
  }
  glLoadMatrixd((double*)mat);
}



void DSRGlWindow::mouseMoveEvent(QMouseEvent *event){
  //! Either rotates (left mouse button) the molecule
  //! or zooms the molecule (right mouse button).
  double x = event->pos().x();
  double y = event->pos().y();
  GLfloat dx = GLfloat( event->x() - moux) / width();
  GLfloat dy = GLfloat( event->y() - mouy) / height();
  if((event->buttons() & Qt::RightButton)){
    glScaled(1.0-dy,1.0-dy,1.0-dy);
  }else{
    glRotateL(dy*360.0,1.0,0.0,0.0);
    glRotateL(dx*360.0,0.0,1.0,0.0);
  }
  moux=x;
  mouy=y;
  updateGL();
} 

void DSRGlWindow::wheelEvent(QWheelEvent *event){
  /*
modifiers:
Qt::NoModifier	        0x00000000  No modifier key is pressed.
Qt::ShiftModifier	      0x02000000  A Shift key on the keyboard is pressed.
Qt::ControlModifier	    0x04000000  A Ctrl key on the keyboard is pressed.
Qt::AltModifier	        0x08000000  An Alt key on the keyboard is pressed.
Qt::MetaModifier	      0x10000000  A Meta key on the keyboard is pressed.
Qt::KeypadModifier	    0x20000000  A keypad button is pressed.
Qt::GroupSwitchModifier	0x40000000  X11 only. A Mode_switch key on the keyboard is pressed.
*/ 
  int numDegrees = event->delta() / 8;
  int numSteps = numDegrees / 15;
  if (event->modifiers()==Qt::NoModifier){
    int d = myFont.pointSize();
    d = (d+numSteps>4)?d+numSteps:d;
    d = qMax(d,7);
    myFont.setPointSize(d);
    updateGL();
  }

}

void DSRGlWindow::tryfit(){
  if (showFit==NULL) return;
  fita.clear();
  fitabonds.clear();
  if ((selFragAt.size()!=3)||(m_molecule->selectedatoms.size()!=3)) {updateGL();return;}
  
  //tree orthogonal vectors of unit size form a rotation matrix 
  //the transponse of the latter is its inverse
  V3 a1,a2,a3;
  a1=Normalize(m_molecule->selectedatoms.at(0).pos-m_molecule->selectedatoms.at(1).pos);//a1=(A-B)/|A-B|
  a2=Normalize(m_molecule->selectedatoms.at(0).pos-m_molecule->selectedatoms.at(2).pos);//a2=(A-C)/|A-B|
  a3=Normalize(a1%a2);//a3=(a1 x a2)/(a1 x a2)
  a2=a3%a1;//a2 = a3 x a1  
  Matrix amat=Matrix(a1,a2,a3);

  V3 b1,b2,b3;
  b1=Normalize(selFragAt.at(0).pos-selFragAt.at(1).pos);
  b2=Normalize(selFragAt.at(0).pos-selFragAt.at(2).pos);
  b3=Normalize(b1%b2);
  b2=b3%b1;
  Matrix bmat=transponse(Matrix(b1,b2,b3));
  if ( (a1*(a2%a3)<0.9) || (b1*(b2%b3)<0.9) ) {  //linear dependence id determinant is 0
    updateGL();
    return;
  }
  // m_molecule->showatoms was m_molecule->asymm before, but showatoms is needed for disorder
  // around special positions:
  for (int i=0; i<m_molecule->showatoms.size(); i++) {
    fita.append(m_molecule->showatoms[i]);
    fita.last().pos+=-1.0*(m_molecule->selectedatoms.at(0).pos); 
    fita.last().pos=fita.last().pos*amat;
    fita.last().pos=fita.last().pos*bmat;
    fita.last().pos+=selFragAt.at(0).pos;
    fita.last().part=1;
    if (Norm (fita.last().pos) > 12) {  // 4 Angstroem cut-off (3*4=12)
      fita.removeLast();
    }
  }  
  fitabonds=mole.connecting(fita,true);
  updateGL();
  printf("FIT:\n================================================================================\n");
  for (int i=0; i<xd.size(); i++){
    int imin=0;
    double dmin=10.0;
    for (int j=0; j<m_molecule->showatoms.size(); j++){
      V3 ort=m_molecule->showatoms.at(j).pos;
      ort+=-1.0*(m_molecule->selectedatoms.at(0).pos); 
      ort=ort*amat;
      ort=ort*bmat;
      ort+=selFragAt.at(0).pos;

      double d=Distance(ort,xd.at(i).pos);
      if (d<dmin) imin=j;
      dmin=(d<dmin)?d:dmin;
    }
    printf("%-8s <==> %-8s %f\n", xd.at(i).Label.toStdString().c_str(),
           m_molecule->showatoms.at(imin).Label.toStdString().c_str(), sqrt(dmin));
  }
}