/*******************************************************************************
*                                                                              *
*                                   Viewmol                                    *
*                                                                              *
*                               M K C O N N . C                                *
*                                                                              *
*                 Copyright (c) Joerg-R. Hill, October 2003                    *
*                                                                              *
********************************************************************************
*
* $Id: mkconn.c,v 1.7 2004/08/29 14:55:28 jrh Exp $
* $Log: mkconn.c,v $
* Revision 1.7  2004/08/29 14:55:28  jrh
* Release 2.4.1
*
* Revision 1.6  2003/11/07 11:07:02  jrh
* Release 2.4
*
* Revision 1.5  2000/12/10 15:11:48  jrh
* Release 2.3
*
* Revision 1.4  1999/05/24 01:26:32  jrh
* Release 2.2.1
*
* Revision 1.3  1999/02/07 21:52:40  jrh
* Release 2.2
*
* Revision 1.2  1998/01/26 00:48:40  jrh
* Release 2.1
*
* Revision 1.1  1996/12/10  18:42:09  jrh
* Initial revision
*
*/
#include<math.h>
#include<stdio.h>
#include<stdlib.h>
#include<X11/Intrinsic.h>
#include<Xm/Xm.h>
#include "viewmol.h"
#include "dialog.h"

extern struct MOLECULE *molecules;
extern struct WINDOW windows[];
extern double amplitude, forceScale, hbondThreshold, lineWidth;
extern int nmolecule, ne;
extern int debug, bondType;
extern Widget topShell;

double centerMolecule(struct MOLECULE *);
void setWindowSize(double);

extern struct SAVE *saveGeometry(struct MOLECULE *);
extern double bondLength(struct ATOM *, int, int);
extern void GetMessageBoxButton(Widget, XtPointer, caddr_t);
extern char *getStringResource(Widget, char *);
extern void *getmem(size_t, size_t);
extern void *expmem(void *, size_t, size_t);
extern void fremem(void **);
extern void bondOrder(struct MOLECULE *, struct BOND *, int *, int);
extern void findConjugation(struct MOLECULE *, struct BOND *, int *, int, int);
extern void calculateShift(struct MOLECULE *, struct BOND *, int *, int, int,
				   int);
extern void pixelToWorld(int, double *, double *);
extern void redraw(int);
extern clock_t getCPUTime(void);

int makeConnectivity(struct MOLECULE *mol, int setWindow, int saveCoord)
{
  Dimension width, height;
  double r1, r2, box, xpix, ypix;
  clock_t startTime=0;
  int *connectivity, maxBondsPerAtom=8, found, nbonds;
  size_t mnb=20;
  register int i, j, k, l;

  box=centerMolecule(mol);
  if (setWindow) setWindowSize(box);
  if (saveCoord)
  {
    if (mol->coord != NULL) fremem((void **)&(mol->coord));
    mol->coord=saveGeometry(mol);
  }

  if (debug == 1) startTime=getCPUTime();
  repeat:
  connectivity=(int *)getmem((size_t)(maxBondsPerAtom*mol->na), sizeof(int));
  for (i=0; i<mol->na; i++)
  {
    mol->atoms[i].nbonds=0;
    mol->atoms[i].hbonds=0;
    for (j=0; j<mol->na; j++)
    {
      if (i != j)
      {
        r1=bondLength(mol->atoms, i, j);
        r2=mol->atoms[i].radScale*mol->atoms[i].rad
          +mol->atoms[j].radScale*mol->atoms[j].rad;
        if (r1 < r2 && strncmp(mol->atoms[i].element->symbol, "??", 2) &&
            strncmp(mol->atoms[j].element->symbol, "??", 2))
        {
          mol->atoms[i].nbonds++;
          if (mol->atoms[i].nbonds > maxBondsPerAtom)
          {
            maxBondsPerAtom+=5;
            fremem((void **)&connectivity);
            goto repeat;
          }
          connectivity[maxBondsPerAtom*i+mol->atoms[i].nbonds-1]=j;
        }
        else
        {
          if (((!strcmp(mol->atoms[i].name, "H") &&
                 strcmp(mol->atoms[j].name, "H") &&
                 strcmp(mol->atoms[j].name, "C")) ||
               (!strcmp(mol->atoms[j].name, "H") &&
                 strcmp(mol->atoms[i].name, "H") &&
                 strcmp(mol->atoms[i].name, "C"))) && r1 < hbondThreshold &&
                 strncmp(mol->atoms[i].element->symbol, "??", 2) &&
                 strncmp(mol->atoms[j].element->symbol, "??", 2))
          {
            mol->atoms[i].nbonds++;
            mol->atoms[i].hbonds++;
            if (mol->atoms[i].nbonds > maxBondsPerAtom)
            {
              maxBondsPerAtom+=5;
              fremem((void **)&connectivity);
              goto repeat;
            }
            connectivity[maxBondsPerAtom*i+mol->atoms[i].nbonds-1]=(-j);
          }
        }
      }
    } 
  }
  if (debug == 1)
  {
    double clockTick;
#if defined LINUX || defined DARWIN
    clockTick=(double)sysconf(_SC_CLK_TCK);
#else
    clockTick=(double)CLK_TCK;
#endif
    printf("Time to find bonds: %f ms\n", 1000.*(getCPUTime()-startTime)/clockTick);
  }
  XtVaGetValues(windows[VIEWER].widget, XtNwidth, &width, XtNheight, &height, NULL);
  if (mol->bondShift == 0.0)
  {
    pixelToWorld(VIEWER, &xpix, &ypix);
    if (lineWidth == 0.0)
      mol->bondShift=0.005*(double)(width+height)*xpix;
    else
      mol->bondShift=lineWidth*xpix;
  }

/* Fill the bonds structure, but check whether user has deleted a particular
   bond. If an atom is only bond to one other atom assign bond order
   immediately since there can be no dependance on any other atom. */

  if (mol->bonds) fremem((void **)&(mol->bonds));
  mol->bonds=(struct BOND *)getmem(mnb, sizeof(struct BOND));
  mol->nb=0;
  for (i=0; i<mol->na; i++)
  {
    nbonds=mol->atoms[i].nbonds;
    for (j=0; j<nbonds; j++)
    {
      k=connectivity[maxBondsPerAtom*i+j];
      if (i < abs(k))
      {
        found=FALSE;
        for (l=0; l<mol->nbDeleted; l++)
        {
          if ((mol->deletedBonds[2*l] == i && mol->deletedBonds[2*l+1] == abs(k)) ||
              (mol->deletedBonds[2*l+1] == i && mol->deletedBonds[2*l] == abs(k)))
          {
            found=TRUE;
            break;
          }
        }
        if (!found)
        {
          mol->bonds[mol->nb].first=i;
          mol->bonds[mol->nb].second=abs(k);
          if (bondType != SINGLE_BONDS)
          {
            if (mol->atoms[i].nbonds == 1 || mol->atoms[abs(k)].nbonds == 1)
              bondOrder(mol, &(mol->bonds[mol->nb]), connectivity,
                        maxBondsPerAtom);
            else
              mol->bonds[mol->nb].order=0;
          }
          else
            mol->bonds[mol->nb].order=1;
          if (k < 0) mol->bonds[mol->nb].order=(-1);
	  /* Use the fraction of the bond length corresponding to the Van der Waals
	     radius of the atom to color a bond according to the atom */
          mol->bonds[mol->nb].frac=mol->atoms[i].rad/(mol->atoms[i].rad+mol->atoms[abs(k)].rad);
          if (++mol->nb >= mnb)
          {
            mnb+=20;
            mol->bonds=(struct BOND *)expmem((void *)mol->bonds, mnb,
                       sizeof(struct BOND));
          }
        }
      }
    }
  }

/* Correct nbonds if bonds have been deleted. This cannot be done before since
   nbonds is used as indicator for the number of entries in connectivity. */

  for (i=0; i<mol->nbDeleted; i++)
  {
    mol->atoms[mol->deletedBonds[2*i]].nbonds--;
    mol->atoms[mol->deletedBonds[2*i+1]].nbonds--; 
  }

/* Add user defined bonds to the bonds array, but check whether that bond
   was already assigned by the previous loop. */

  for (i=0; i<mol->nbAdded; i++)
  {
    found=FALSE;
    for (j=0; j<mol->nb; j++)
    {
      if ((mol->addedBonds[3*i] == mol->bonds[j].first &&
           mol->addedBonds[3*i+1] == mol->bonds[j].second) ||
          (mol->addedBonds[3*i] == mol->bonds[j].second &&
           mol->addedBonds[3*i+1] == mol->bonds[j].first))
      {
        found=TRUE;
        mol->bonds[j].order=mol->addedBonds[3*i+2];
        if (mol->bonds[j].order > 1)
          calculateShift(mol, &mol->bonds[j], connectivity, maxBondsPerAtom,
                         mol->bonds[j].first, mol->bonds[j].second);
        break;
      }
    }
    if (!found)
    {
      mol->atoms[mol->addedBonds[3*i]].nbonds++;
      mol->atoms[mol->addedBonds[3*i+1]].nbonds++;
      mol->bonds[mol->nb].first=mol->addedBonds[3*i];
      mol->bonds[mol->nb].second=mol->addedBonds[3*i+1];
      mol->bonds[mol->nb].order=mol->addedBonds[3*i+2];
      if (mol->addedBonds[3*i+2] > 1)
        calculateShift(mol, &mol->bonds[mol->nb], connectivity, maxBondsPerAtom,
                       mol->addedBonds[3*i], mol->addedBonds[3*i+1]);
      mol->bonds[mol->nb++].frac=mol->atoms[mol->addedBonds[3*i]].rad
                                /bondLength(mol->atoms, mol->addedBonds[3*i],
                                mol->addedBonds[3*i+1]);
      if (mol->nb >= mnb)
      {
        mnb+=20;
        mol->bonds=(struct BOND *)expmem((void *)mol->bonds, mnb,
                                         sizeof(struct BOND));
      }
    }
  }

/* Assign bond orders for the remaining bonds */

  for (i=0; i<mol->nb; i++)
  {
    if (mol->bonds[i].order == 0)
    {
      if (bondType != SINGLE_BONDS)
        bondOrder(mol, &mol->bonds[i], connectivity, maxBondsPerAtom);
      else
        mol->bonds[i].order=1;
    }
  }

/* Check if there are conjugated double bonds */

  if (bondType == CONJUGATION) findConjugation(mol, mol->bonds,
                                 connectivity, maxBondsPerAtom, mol->nb);

  fremem((void **)&connectivity);
  if (mol->nb > 0)
    mol->bonds=(struct BOND *)expmem((void *)mol->bonds,
                          mol->nb, sizeof(struct BOND));

/* Add unit cell edges */

  if (mol->unitcell)
  {
    mol->unitcell->edges=(struct EDGE *)getmem(12, sizeof(struct EDGE));
    mol->unitcell->ne=0;
    mol->unitcell->edges[mol->unitcell->ne].first=0;
    mol->unitcell->edges[mol->unitcell->ne++].second=1;
    mol->unitcell->edges[mol->unitcell->ne].first=0;
    mol->unitcell->edges[mol->unitcell->ne++].second=2;
    mol->unitcell->edges[mol->unitcell->ne].first=0;
    mol->unitcell->edges[mol->unitcell->ne++].second=3;
    mol->unitcell->edges[mol->unitcell->ne].first=1;
    mol->unitcell->edges[mol->unitcell->ne++].second=4;
    mol->unitcell->edges[mol->unitcell->ne].first=1;
    mol->unitcell->edges[mol->unitcell->ne++].second=7;
    mol->unitcell->edges[mol->unitcell->ne].first=2;
    mol->unitcell->edges[mol->unitcell->ne++].second=6;
    mol->unitcell->edges[mol->unitcell->ne].first=2;
    mol->unitcell->edges[mol->unitcell->ne++].second=7;
    mol->unitcell->edges[mol->unitcell->ne].first=3;
    mol->unitcell->edges[mol->unitcell->ne++].second=4;
    mol->unitcell->edges[mol->unitcell->ne].first=3;
    mol->unitcell->edges[mol->unitcell->ne++].second=6;
    mol->unitcell->edges[mol->unitcell->ne].first=4;
    mol->unitcell->edges[mol->unitcell->ne++].second=5;
    mol->unitcell->edges[mol->unitcell->ne].first=5;
    mol->unitcell->edges[mol->unitcell->ne++].second=6;
    mol->unitcell->edges[mol->unitcell->ne].first=5;
    mol->unitcell->edges[mol->unitcell->ne++].second=7;
  }
  if (debug == 1)
  {
    printf("Bonds:\n");
    for (i=0; i<mol->nb; i++)
    {
      printf("%d: %s(%d)-%s(%d) = %d, frac = %f", i, mol->atoms[mol->bonds[i].first].name,
             mol->bonds[i].first+1, mol->atoms[mol->bonds[i].second].name,
             mol->bonds[i].second+1, mol->bonds[i].order, mol->bonds[i].frac);
      if (abs(mol->bonds[i].order) != 1)
        printf(", multiple bond shift %f %f %f\n", mol->bonds[i].x,
               mol->bonds[i].y, mol->bonds[i].z);
      else
        printf("\n");
    }
    if (mol->unitcell)
    {
      printf("Unit cell edges:\n");
      for (i=0; i<mol->unitcell->ne; i++)
        printf("%d: %d-%d\n", i, mol->bonds[i].first+1, mol->bonds[i].second+1);
    }
  }
  return(TRUE);
}

double centerMolecule(struct MOLECULE *mol)
{
  double mass, box, r;
  int i;

  /* Translate molecule so that it is centered at its
     center of gravity and determine the size of the cube
     enclosing it taking into account the van der Waals
     radii of the atoms */
  mol->transx=0.0;
  mol->transy=0.0;
  mol->transz=0.0;
  mass=0.0;
  for (i=0; i<mol->na; i++)
  {
    mol->transx+=mol->atoms[i].mass*mol->atoms[i].x;
    mol->transy+=mol->atoms[i].mass*mol->atoms[i].y;
    mol->transz+=mol->atoms[i].mass*mol->atoms[i].z;
    mass+=mol->atoms[i].mass;
  }
  mol->transx=-mol->transx/mass;
  mol->transy=-mol->transy/mass;
  mol->transz=-mol->transz/mass;
  if (debug == 1)
    printf("Shifting molecule by %f %f %f.\n", mol->transx, mol->transy,
           mol->transz);
  box=0.0;
  for (i=0; i<mol->na; i++)
  {
    mol->atoms[i].x+=mol->transx;
    mol->atoms[i].y+=mol->transy;
    mol->atoms[i].z+=mol->transz;
    r=mol->atoms[i].radScale*mol->atoms[i].rad;
    box=fabs(mol->atoms[i].x)+r > box ? fabs(mol->atoms[i].x)+r : box;
    box=fabs(mol->atoms[i].y)+r > box ? fabs(mol->atoms[i].y)+r : box;
    box=fabs(mol->atoms[i].z)+r > box ? fabs(mol->atoms[i].z)+r : box;
  }
  if (mol->unitcell)
  {
    for (i=0; i<mol->unitcell->nc; i++)
    {
      mol->unitcell->corners[i].x+=mol->transx;
      mol->unitcell->corners[i].y+=mol->transy;
      mol->unitcell->corners[i].z+=mol->transz;
      r=mol->unitcell->corners[i].radScale*mol->unitcell->corners[i].rad;
      box=fabs(mol->unitcell->corners[i].x)+r > box ? fabs(mol->unitcell->corners[i].x)+r : box;
      box=fabs(mol->unitcell->corners[i].y)+r > box ? fabs(mol->unitcell->corners[i].y)+r : box;
      box=fabs(mol->unitcell->corners[i].z)+r > box ? fabs(mol->unitcell->corners[i].z)+r : box;
    }
  }
  if (box == 0.0) box=mol->atoms[0].rad;
  /* Add 5 % space and make it the space diagonal of the cube
     so that rotations do not move any atom out of the clipping
     volume */
  box*=1.05*sqrt(3.);
  return(box);
}

void setBondType(Widget widget, caddr_t which, XmToggleButtonCallbackStruct *data)
{
  register int i;

  if (data->set)
  {
    bondType=(int)which;
    for (i=0; i<nmolecule; i++)
      makeConnectivity(&molecules[i], FALSE, FALSE);
    redraw(VIEWER);
  }
}

void setWindowSize(double box)
{
  Dimension width, height;
  double r;

  amplitude=forceScale=box;
  windows[VIEWER].left=(-box);
  windows[VIEWER].right=box;
  windows[VIEWER].bottom=(-box);
  windows[VIEWER].top=box;
  windows[VIEWER].near=(-box);
  windows[VIEWER].far=box;

  XtVaGetValues(windows[VIEWER].widget, XtNwidth, &width, XtNheight, &height, NULL);
  if (width > height)
  {
    r=((double)width/(double)height-1.0)*windows[VIEWER].near;
    windows[VIEWER].left=windows[VIEWER].near+r;
    windows[VIEWER].right=windows[VIEWER].far-r;
    windows[VIEWER].bottom=windows[VIEWER].near;
    windows[VIEWER].top=windows[VIEWER].far;
  }
  else
  {
    r=((double)height/(double)width-1.0)*windows[VIEWER].near;             
    windows[VIEWER].left=windows[VIEWER].near;                               
    windows[VIEWER].right=windows[VIEWER].far;
    windows[VIEWER].bottom=windows[VIEWER].near+r;                         
    windows[VIEWER].top=windows[VIEWER].far-r;
  }
}
/*
void deleteBond(struct MOLECULE *mol, int i, int j)
{
  register int i;

  for (i=0; i<mol->nb; i++)
  {
    if ((molecules[imol].bonds[i].first == clicked[0] && molecules[imol].bonds[i].second == clicked[1])
      ||(molecules[imol].bonds[i].first == clicked[1] && molecules[imol].bonds[i].second == clicked[0]))
    {
      for (j=i+1; j<molecules[imol].nb; j++)
      {
                                      molecules[imol].bonds[j-1].first=molecules[imol].bonds[j].first;
                                      molecules[imol].bonds[j-1].second=molecules[imol].bonds[j].second;
                                      molecules[imol].bonds[j-1].order=molecules[imol].bonds[j].order;
                                      molecules[imol].bonds[j-1].frac=molecules[imol].bonds[j].frac;
                                      molecules[imol].bonds[j-1].x=molecules[imol].bonds[j].x;
                                      molecules[imol].bonds[j-1].y=molecules[imol].bonds[j].y;
                                      molecules[imol].bonds[j-1].z=molecules[imol].bonds[j].z;
                                    }
                                    molecules[imol].nb--;
                                    molecules[imol].bonds=(struct BOND *)expmem((void *)molecules[imol].bonds,
                                                            molecules[imol].nb, sizeof(struct BOND));
                                  } */