RigidCLL.cpp

Go to the documentation of this file.

#include "RigidCLL.h"

#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <limits>

#define BB_MARGIN 2

#define MAX_UINT std::numeric_limits<unsigned int>::max()

#define ERROR(s) {cerr<<s<<"\n";exit(0);}

#define INF std::numeric_limits<double>::max()

#define MANHATTAN_DISTANCE(dx,dy,dz,a1,a2)  dx=a1->cp[0]-a2->cp[0];dy=a1->cp[1]-a2->cp[1];dz=a1->cp[2]-a2->cp[2]

/*************************************************************************/
/*************************************************************************/
/*                          AtomInfo methods                             */
/*************************************************************************/
/*************************************************************************/

AtomInfo::AtomInfo(unsigned int identifier,unsigned int rigid,double x,double y,double z,double radius)
{
  if (radius<=0)
    ERROR("The atom radius must be positive.");
  
  id=identifier;
  rID=rigid;

  /* The current and original positions of the atom */
  cp[0]=op[0]=x;
  cp[1]=op[1]=y;
  cp[2]=op[2]=z;

  r=radius;

  /* Indexes of the cell including this atom */
  c[0]=c[1]=c[2]=0; /* not used until the grid is defined. */
}

void AtomInfo::move(RigidTransform t)
{
  cp[0]=t[0][0]*op[0]+t[0][1]*op[1]+t[0][2]*op[2]+t[0][3];
  cp[1]=t[1][0]*op[0]+t[1][1]*op[1]+t[1][2]*op[2]+t[1][3];
  cp[2]=t[2][0]*op[0]+t[2][1]*op[1]+t[2][2]*op[2]+t[2][3];
}

double AtomInfo::distance(AtomInfo *a)
{
  double dx,dy,dz;

  MANHATTAN_DISTANCE(dx,dy,dz,this,a);

  return(sqrt(dx*dx+dy*dy+dz*dz));
}

ostream& operator<<(ostream& os,const AtomInfo &a)
{
  os << "id:" << a.id ;
  os << " r:" << a.rID ;
  os << " op:(" << a.op[0] << "," << a.op[1] << "," << a.op[2] << ")";
  os << " cp:(" << a.cp[0] << "," << a.cp[1] << "," << a.cp[2] << ")";

  return(os);
}

/*************************************************************************/
/*************************************************************************/
/*                          AtomPair methods                             */
/*************************************************************************/
/*************************************************************************/

AtomPair::AtomPair(AtomInfo *atom1,AtomInfo *atom2,double dst2)
{
  a1=atom1->id;
  r1=atom1->rID;
  a2=atom2->id;
  r2=atom2->rID;
  d2=dst2;
}

bool AtomPair::operator<(const AtomPair a) const
{
  unsigned int min1,min2,max1,max2;

  min1=(this->a1<this->a2?this->a1:this->a2); // minimum of the two atoms in p1
  max1=(this->a1<this->a2?this->a2:this->a1); // maximum of the two atoms in p1

  min2=(a.a1<a.a2?a.a1:a.a2); // minimum of the two atoms in p2
  max2=(a.a1<a.a2?a.a2:a.a1); // maximum of the two atoms in p2

  return((min1<min2)||((min1==min2)&&(max1<max2)));
}

ostream& operator<<(ostream& os,const AtomPair &a)
{
  if (a.a1<a.a2)
    os << "  (" << a.a1 << "[" << a.r1  << "]," << a.a2 << "[" << a.r2  << "])";  
  else
    os << "  (" << a.a2 << "[" << a.r2  << "]," << a.a1 << "[" << a.r1  << "])";

  os << " d: " << sqrt(a.d2);  
  return(os);
}

bool AtomPair::operator==(AtomPair &a)
{
  /* we only need to compare atom identifiers */
  return(((this->a1==a.a1)&&(this->a2==a.a2))||
         ((this->a1==a.a2)&&(this->a2==a.a1)));
}

/*************************************************************************/
/*************************************************************************/
/*                          CellInfo methods                             */
/*************************************************************************/
/*************************************************************************/

CellInfo::CellInfo()
{
  mobileID=MAX_UINT;
  fixedID=MAX_UINT;
}

void CellInfo::Reset()
{
  mobileID=MAX_UINT;
}

void CellInfo::AddFixedAtom(unsigned int id,unsigned int *fixedLL)
{
  fixedLL[id]=fixedID;
  fixedID=id;
}

void CellInfo::AddMobileAtom(unsigned int id,unsigned int *mobileLL)
{
  mobileLL[id]=mobileID;
  mobileID=id;
}

/*************************************************************************/
/*************************************************************************/
/*                     RigidCLL private methods                          */
/*************************************************************************/
/*************************************************************************/

void RigidCLL::InitRigidCLL()
{
  o[0]=o[1]=o[2]=0.0;
  s[0]=s[1]=s[2]=0;
  grid=NULL;
  fixedLL=NULL;
  mobileLL=NULL;
  cutoff=0.0; /* default is steric clash */
  d=0.0; /* updated in warmup */
  cr.firstAtom=cr.lastAtom=MAX_UINT;
  nFixed=0;
  nMobile=0;
  rID=0;
  lastInGrid=MAX_UINT;
  checkInteraction=NULL;
}

void RigidCLL::UpdateCell(AtomInfo *a)
{
  unsigned int i;
  
  for(i=0;i<3;i++)
    {
      a->c[i]=(unsigned int)floor((a->cp[i]-o[i])/d);
      #if (RIGIDCLL_DEBUG)
        /* atoms can not be in the limit of the grid, otherwise
           the search will fail */
        if ((a->c[i]<=0)||(a->c[i]>=s[i]-1))
          ERROR("Cell out of range in UpdateCell");
      #endif
    }
}

void RigidCLL::InsertFixedInGrid(unsigned int n)
{
  AtomInfo *a;
  
  #if (RIGIDCLL_DEBUG)
    if (grid==NULL)
      ERROR("Inserting atom in an empty grid. Use the WarmUp method.");
  #endif

  a=&(fixedAtoms[n]);

  /* For fixed atoms we do not have the grid cell pre-computed */
  UpdateCell(a);

  grid[a->c[0]][a->c[1]][a->c[2]].AddFixedAtom(n,fixedLL);

  /* fixed atoms are never removed from the grid: there is no need
     to update the 'lastInGrid' */
}

void RigidCLL::InsertMobileInGrid(unsigned int n)
{
  AtomInfo *a;
  
  #if (RIGIDCLL_DEBUG)
    if (grid==NULL)
      ERROR("Inserting atom in an empty grid. Use the WarmUp method.");
  #endif

  /* for mobile atoms the position in the grid is already pre-computed */

  a=&(mobileAtoms[n]);

  grid[a->c[0]][a->c[1]][a->c[2]].AddMobileAtom(n,mobileLL);
  
  /* For mobile atoms we have to keep track of the last
     atom added to the grid. We assume that 'n' is
     monotonically increasing. */
  lastInGrid=n;
}

void RigidCLL::DeleteGrid()
{
  if (grid!=NULL)
    {
      /* It would be a good idea to use the atoms (fixed and mobile) to identify
         the possibly used cells in the grid. In any case, to deallocate we
         need the two loos 'i' and 'j'.... so adding a fast loop over 'k'
         is probably a reasonable solution. */

      unsigned int i,j;

      for(i=0;i<s[0];i++)
        {
          for(j=0;j<s[1];j++)
            {
              delete[] grid[i][j];
            }
          delete[] grid[i];
        }
      delete[] grid;

      delete[] fixedLL;
      delete[] mobileLL;
      
      grid=NULL;
      fixedLL=NULL;
      mobileLL=NULL;
    }
}

void RigidCLL::DeleteInteractions()
{
  if (checkInteraction!=NULL)
    {
      unsigned int i,nr;

      nr=rID+1;
      for(i=0;i<nr;i++)
        {
          delete [] checkInteraction[i];
        }
      delete [] checkInteraction;

      checkInteraction=NULL;
    }
}

void RigidCLL::InitInteractions()
{
  unsigned int i,j,nr;

  /* If defining a rigid, close it first */
  if (cr.firstAtom!=MAX_UINT)
    EndRigid();

  /* now define the interaction table */
  nr=rID+1;
  checkInteraction=new bool*[nr];
  for(i=0;i<nr;i++)
    {
      checkInteraction[i]=new bool[nr];
      for(j=0;j<nr;j++)
        checkInteraction[i][j]=(i!=j); /* self-interactions are not considered */
    }
}


void RigidCLL::ResetQuery()
{
  unsigned int j;
  AtomInfo *a1;

  if (lastInGrid!=MAX_UINT)
    {
      for(j=0;j<=lastInGrid;j++)
        {
          a1=&(mobileAtoms[j]);
          grid[a1->c[0]][a1->c[1]][a1->c[2]].Reset();
        }    
      lastInGrid=MAX_UINT;
    }
}

void RigidCLL::Move(RigidTransform *t)
{
  unsigned int j;
  AtomInfo *a1;
      
  /* Move the atoms */
  for(j=0;j<nMobile;j++)
    {
      a1=&(mobileAtoms[j]);
      a1->move(t[a1->rID]);
    }
}

bool RigidCLL::Check27(AtomInfo *a1)
{
  AtomInfo *a2;
  unsigned int ix,iy,iz,aID;
  double dx,dy,dz,r;
  bool collision=false;

  /* for neighbours */
  ix=a1->c[0]-1;
  while((!collision)&&(ix<=a1->c[0]+1)) 
    {
      iy=a1->c[1]-1;
      while((!collision)&&(iy<=a1->c[1]+1))
        {              
          iz=a1->c[2]-1;
          while((!collision)&&(iz<=a1->c[2]+1))
            {
              /* first we check for possible interaction with
                 the fixed atoms in the cell, if any */
              if (checkInteraction[a1->rID][0])
                {
                  /* loop over the list of fixed atoms */
                  aID=grid[ix][iy][iz].fixedID;
                  while((!collision)&&(aID!=MAX_UINT))
                    {
                      a2=&(fixedAtoms[aID]);
                                
                      MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                      r=a1->r+a2->r;

                      collision=(dx*dx+dy*dy+dz*dz<=r*r);
                              
                      aID=fixedLL[aID];
                    }
                }
              /* and now loop over the mobile atoms, if any */
              aID=grid[ix][iy][iz].mobileID;
              while((!collision)&&(aID!=MAX_UINT))
                {
                  a2=&(mobileAtoms[aID]);
                                              
                  if (checkInteraction[a1->rID][a2->rID]) 
                    {
                      MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                      r=a1->r+a2->r;

                      collision=(dx*dx+dy*dy+dz*dz<=r*r);
                    }
                  aID=mobileLL[aID];
                }
              iz++;
            }
          iy++;
        }
      ix++;
    }
  
  return(collision);
}

/*************************************************************************/
/*************************************************************************/
/*                      RigidCLL public methods                          */
/*************************************************************************/
/*************************************************************************/

RigidCLL::RigidCLL():fixedAtoms(),mobileAtoms(),rigids()
{
  InitRigidCLL();
}

RigidCLL::RigidCLL(char *name):fixedAtoms(),mobileAtoms(),rigids()
{
  ifstream file;
  unsigned int nr,na,i,j,k;
  double x,y,z,r;
  bool f;

  /* First, the standard init */
  InitRigidCLL();

  /* Now, we fill in the structure using data from the file */
  file.open(name);

  file>>nr;

  k=0;
  for(i=0;i<nr;i++)
    {
      file >> na;
      file >> f;

      if (!f)
        NewRigid();
      
      for(j=0;j<na;j++)
        {
          file >> x;
          file >> y;
          file >> z;
          file >> r;

          if (f)
            AddFixedAtom(k,x,y,z,r);
          else
            AddAtom(k,x,y,z,r);
          k++;
        }

      if (!f)
        EndRigid();
    }

  file.close();

  /* allow the interaction between all the rigid groups */
  InitInteractions();
}

unsigned int RigidCLL::GetNumRigids()
{
  /* The group of fixed atoms counts as one rigid */
  return(rID+1);
}

void RigidCLL::SetCutOff(double c)
{
  if (c<0.0)
    ERROR("The cutoff must be possitive");

  cutoff=c;

  /* a new cutoff -> re-define the grid. We delete
     it and it will be re-defined when necessary. */
  DeleteGrid();
}

double RigidCLL::GetCutOff()
{  
  return(cutoff);
}

double RigidCLL::GetBBSize(unsigned int axis)
{
  if ((axis<0)||(axis>2))
    ERROR("Axis out of range in GetBBSize");

  if ((grid==NULL)||(s[axis]==0))
    ERROR("Calling GetBBSize on non-initialized RigidCLL object (call the WarmUp method first)");

  return(s[axis]*d);
}

RigidCLL::~RigidCLL()
{
  Reset();
}

void RigidCLL::AddFixedAtom(unsigned int id,double x,double y,double z,double r)
{
  if (checkInteraction!=NULL)
    ERROR("Adding atoms to a closed model. Use the Reset method first.");

  fixedAtoms.push_back(AtomInfo(id,0,x,y,z,r));
  nFixed++;
}

void RigidCLL::AddFixedAtom(unsigned int id,double *p,double r)
{
  AddFixedAtom(id,p[0],p[1],p[2],r);
}

unsigned int RigidCLL::NewRigid()
{
  if (checkInteraction!=NULL)
    ERROR("Adding atoms to a closed model. Use the Reset method first.");

  /* If already defining a rigid, close it first */
  if (cr.firstAtom!=MAX_UINT)
    EndRigid();

  /* initialize the current rigid structure */
  cr.firstAtom=nMobile;
  rID++; /* mobile rigids are numbered form 1 */

  return(rID);
}

void RigidCLL::AddAtom(unsigned int id,double x,double y,double z,double r)
{
  if (cr.firstAtom==MAX_UINT)
    ERROR("Adding an atom out of a rigid group. Use NewRigid first.");

  mobileAtoms.push_back(AtomInfo(id,rID,x,y,z,r));
  nMobile++;
}

void RigidCLL::AddAtom(unsigned int id,double *p,double r)
{
  AddAtom(id,p[0],p[1],p[2],r);
}

void RigidCLL::EndRigid()
{
  /* if we are not defining a rigid, just ignore this call */
  if (cr.firstAtom!=MAX_UINT)
    {
      if (checkInteraction!=NULL)
        ERROR("Modifying a closed model. Use the Reset method first.");
      
      cr.lastAtom=nMobile; /* this is actually the first atom of the next rigid, if any */
      rigids.push_back(cr);
      cr.firstAtom=cr.lastAtom=MAX_UINT; /* flag to indicate the end of the rigid definition */
    }
}

void RigidCLL::ActivateInteraction(unsigned int r1,unsigned int r2)
{
  if (checkInteraction==NULL)
    InitInteractions();

  /* rID is the identifier of the last rigid in the system */
  if ((r1>rID)||(r2>rID))
    ERROR("Rigid cluster identifier out of range in ActivateInteraction");

  if (r1==r2)
    ERROR("Can not activate the self-interaction of a rigid group.");

  checkInteraction[r1][r2]=checkInteraction[r2][r1]=true;
}

void RigidCLL::DeactivateInteraction(unsigned int r1,unsigned int r2)
{
  if (checkInteraction==NULL)
    InitInteractions();

  if ((r1>rID)||(r2>rID))
    ERROR("Rigid cluster identifier out of range in DeactivateInteraction");

  checkInteraction[r1][r2]=checkInteraction[r2][r1]=false;
}

void RigidCLL::WarmUp(RigidTransform *t)
{
  unsigned int i,j;
  double minC[3],maxC[3],maxR,range;
  vector<AtomInfo>::iterator it;

  if (t==NULL)
    ERROR("The Warmup methods need a valid conformation");
  
  /* if the interactions are not defined, assume all rigids interact between them */
  if (checkInteraction==NULL)
    InitInteractions();

  /* If we already have a grid, delete it before re-initializing */
  DeleteGrid();
  
  /* Compute the bounding box for all atoms */
  minC[0]=minC[1]=minC[2]=INF;
  maxC[0]=maxC[1]=maxC[2]=-INF;
  maxR=0;

  /* BB for the fixed atoms (for fixed atoms 'op' = 'cp' always) */
  for(it=fixedAtoms.begin();it!=fixedAtoms.end();++it)
    {
      for(j=0;j<3;j++)
        {
          if (it->cp[j]>maxC[j]) maxC[j]=it->cp[j];
          else { if (it->cp[j]<minC[j]) minC[j]=it->cp[j]; }
        }

      if (it->r>maxR) maxR=it->r;
    }

  /* include the mobile atoms in the BB. Use the current position 'cp'
     and not the origina one 'op' which is in local coordinates many times */
  for(it=mobileAtoms.begin();it!=mobileAtoms.end();++it)
    {
      /* update the 'cp' */
      it->move(t[it->rID]);
      for(j=0;j<3;j++)
        {
          if (it->cp[j]>maxC[j]) maxC[j]=it->cp[j];
          else { if (it->cp[j]<minC[j]) minC[j]=it->cp[j]; }
        }
      if (it->r>maxR) maxR=it->r;
    }

  /* compute the cutoff distance, if not given */
  if (cutoff==0)
    d=2.0*maxR;
  else
    d=cutoff;
  d2=d*d;

  /* Set up the origin of the grid and its size (in number of cells) 
     We define a grid that covers 27 times the bounding box previously
     computed (the bounding box is all surrounded by other boxes of
     the same size). In this way we hope that the motions of the molecule
     would neve scape the grid.
   */
  for(j=0;j<3;++j)
    {
      range=maxC[j]-minC[j];
      o[j]=minC[j]-(BB_MARGIN)*range;
      s[j]=(unsigned int)ceil((1+2*BB_MARGIN)*range/d);
    }

  /* Set up the grid */
  grid=new CellInfo **[s[0]];
  for(i=0;i<s[0];i++)
    {
      grid[i]=new CellInfo *[s[1]];
      for(j=0;j<s[1];j++)
        {
          grid[i][j]=new CellInfo [s[2]];
        }
    }
  fixedLL=new unsigned int[nFixed];
  mobileLL=new unsigned int[nMobile];

  /* Insert the fixed atoms in the grid */
  for(i=0;i<nFixed;i++)
    InsertFixedInGrid(i);

  /* we just created the grid: there is no mobile atom in it yet. */
  lastInGrid=MAX_UINT;
}
              
void RigidCLL::Reset()
{
  DeleteGrid();
  DeleteInteractions();
}

bool RigidCLL::StericClash(RigidTransform *t)
{
  unsigned int i,j;
  AtomInfo *a;
  bool collision=false;

    
  /* We could use the checkInteraction table to more accurately determine wether or not it 
     is worth proceed with the interaction detection. Here we use a conservative condition. */
  if ((rID>1)||((rID==1)&&(nFixed>0)))
    {  
      #if (RIGIDCLL_DEBUG)
        if (cutoff>0)
          ERROR("Using StericClash in a RigidCLL initialized to detect nearby-interactions");
      #endif 

      /* ensure we have a grid and that it is in it's iniital state */
      if (grid==NULL)
        WarmUp(t);
      else
        ResetQuery();

      /* Proceed with the query */
      for(i=0;((!collision)&&(i<rID));i++)
        {
          for(j=rigids[i].firstAtom;((!collision)&&(j<rigids[i].lastAtom));j++)
            {
              a=&(mobileAtoms[j]);

              /* Move the atom (we delay the movement as much as possible) */
              a->move(t[a->rID]); /* if we used the WarmUp this is redudant but... */

              /* This pre-computes the cell for this atom. This must be done even
                 if we skip the query phase (it is used in the insert). */
              UpdateCell(a);

              /* If we do not have fixed atoms, we can skip the query phase for 
                 the first rigid (the grid is empty!). */
              if ((nFixed>0)||(i>0))
                collision=Check27(a);
            }

          /* Insert atoms in the grid (there is no need to add last rigid) */
          if ((!collision)&&(i<rID-1))
            {
              for(j=rigids[i].firstAtom;j<rigids[i].lastAtom;j++)
                InsertMobileInGrid(j);   // this assumes UpdateCell already executed
            }
        }
    }
  return(collision);
}

bool RigidCLL::CLLStericClash(RigidTransform *t)
{
  unsigned int i;
  bool collision=false;

  if ((rID>1)||((rID==1)&&(nFixed>0)))
    {
      #if (RIGIDCLL_DEBUG)
        if (cutoff>0)
          ERROR("Using CLLStericClash in a RigidCLL initialized to detect nearby-interactions");
      #endif
 
      /* ensure we have a grid and that it is in it's iniital state */
      if (grid==NULL)
        WarmUp(t); /* this already moves the atoms */
      else
        {
          ResetQuery();
          Move(t);
        }

      /* Proceed with the query */

      /* First insert all the mobile atoms in the grid */
      for(i=0;i<nMobile;i++)
        {
          UpdateCell(&(mobileAtoms[i]));
          InsertMobileInGrid(i);
        }
  
      /* Now check for collision */
      for(i=0;((!collision)&&(i<nMobile));i++)
        collision=Check27(&(mobileAtoms[i]));
    }

  return(collision);
}

bool RigidCLL::BruteForceStericClash(RigidTransform *t)
{
  unsigned int i,j;
  AtomInfo *a1,*a2;
  double dx,dy,dz,r;
  bool collision=false;

  #if (RIGIDCLL_DEBUG)
    if (cutoff>0)
      ERROR("Using BruteForceStericClash in a RigidCLL initialized to detect nearby-interactions");
  #endif

  if ((rID>1)||((rID==1)&&(nFixed>0)))
    {
      /* ensure the interaction table is defined */
      if (checkInteraction==NULL)
        InitInteractions();
      
      /* Move all atoms */
      Move(t);

      for(i=0;((!collision)&&(i<nMobile));i++)
        {
          a1=&(mobileAtoms[i]);

          if (checkInteraction[a1->rID][0])
            {
              for(j=0;((!collision)&&(j<nFixed));j++)
                {
                  a2=&(fixedAtoms[j]);
                     
                  MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                  r=a1->r+a2->r;
                
                  collision=(dx*dx+dy*dy+dz*dz<=r*r);
                }
            }
        
          for(j=i+1;((!collision)&&(j<nMobile));j++)
            {
              a2=&(mobileAtoms[j]);
                     
              if (checkInteraction[a1->rID][a2->rID])
                {
                  MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                  r=a1->r+a2->r;
              
                  collision=(dx*dx+dy*dy+dz*dz<=r*r);
                }
            }
        }
    }

  return(collision);
}

void RigidCLL::Interactions(RigidTransform *t,list<AtomPair> *c)
{
  unsigned int i,j,aID,ix,iy,iz;
  double dx,dy,dz,d;
  AtomInfo *a1,*a2;
  RigidInfo *rigid;
    
  /* We could use the checkInteraction table to more accurately determine wether or not it 
     is worth proceed with the interaction detection. Here we use a conservative condition. */
  if ((rID>1)||((rID==1)&&(nFixed>0)))
    {
      #if (RIGIDCLL_DEBUG)
        if (cutoff==0)
          ERROR("Using Interactions in a RigidCLL initialized to detect steric clashes");
      #endif

      /* ensure we have a grid and that it is in it's iniital state */
      if (grid==NULL)
        WarmUp(t); /* this already moves the atoms */
      else
        {
          ResetQuery();
          Move(t);
        }

      /* Proceed with the query */
      for(i=0;i<rID;i++)
        {
          rigid=&(rigids[i]);

          for(j=rigid->firstAtom;j<rigid->lastAtom;j++) /* atom in rigid */
            {
              a1=&(mobileAtoms[j]);

              /* This pre-computes the cell for this atom and must be done even
                 if we skip the query phase. */
              UpdateCell(a1); 

              /* If we do not have fixed atoms, we can skip the query phase for 
                 the first rigid (the grid is empty!) */
              if ((nFixed>0)||(i>0))
                {
                  /* for neighbours */
                  for(ix=a1->c[0]-1;ix<=a1->c[0]+1;ix++)
                    {
                      for(iy=a1->c[1]-1;iy<=a1->c[1]+1;iy++)
                        {              
                          for(iz=a1->c[2]-1;iz<=a1->c[2]+1;iz++)
                            {
                              if (checkInteraction[a1->rID][0])
                                {
                                  /* first we check for possible interaction with
                                     the fixed atoms in the cell, if any */
                                  aID=grid[ix][iy][iz].fixedID;
                                  while(aID!=MAX_UINT)
                                    {
                                      a2=&(fixedAtoms[aID]);

                                      MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                                      d=dx*dx+dy*dy+dz*dz;
                           
                                      if (d<=d2)
                                        c->push_back(AtomPair(a1,a2,d));

                                      aID=fixedLL[aID];
                                    }
                                }

                              /* and now with the mobile atoms, if any */
                              aID=grid[ix][iy][iz].mobileID;
                              while(aID!=MAX_UINT)
                                {
                                  a2=&(mobileAtoms[aID]);
                               
                                  if (checkInteraction[a1->rID][a2->rID])
                                    {
                                      MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                                      d=dx*dx+dy*dy+dz*dz;
                                          
                                      if (d<=d2)
                                        c->push_back(AtomPair(a1,a2,d));
                                    }
                                  aID=mobileLL[aID];
                                }
                            }
                        }
                    }
                }
            }
      
          /* There is no need to add last rigid group to the grid (would no be used) */
          if (i<rID-1) 
            {
              /* Insert atoms in the grid */
              for(j=rigid->firstAtom;j<rigid->lastAtom;j++)
                InsertMobileInGrid(j);
            }
        }
    }
}

void RigidCLL::CLLInteractions(RigidTransform *t,list<AtomPair> *c)
{
  unsigned int i,aID,ix,iy,iz;
  double dx,dy,dz,d;
  AtomInfo *a1,*a2;

  /* If we have more than one atom -> we have to check for interaction. */
  if ((nMobile>1)||((nMobile==1)&&(nFixed>0)))
    {
      #if (RIGIDCLL_DEBUG)
        if (cutoff==0)
          ERROR("Using CLLInteractions in a RigidCLL initialized to detect steric clashes");
      #endif

      /* ensure we have a grid and that it is in it's iniital state */
      if (grid==NULL)
        WarmUp(t); /* this already moves the atoms */
      else
        {
          ResetQuery();
          Move(t);
        }

      /* Insert all atoms in the grid */
      for(i=0;i<nMobile;i++)
        {
          UpdateCell(&(mobileAtoms[i])); 
          InsertMobileInGrid(i);
        }

      /* Check for all possible interactions */
       for(i=0;i<nMobile;i++)
         {
           a1=&(mobileAtoms[i]);

           /* The loops below are the same as those in MoveAndDetect
               with a minor difference (see comments below).
               This minor difference avoids re-using the same code 
               (as we did in MoveAndCheck)*/
           for(ix=a1->c[0]-1;ix<=a1->c[0]+1;ix++)
             {
               for(iy=a1->c[1]-1;iy<=a1->c[1]+1;iy++)
                 {             
                   for(iz=a1->c[2]-1;iz<=a1->c[2]+1;iz++)
                     {
                       /* DIFFERENCE: do not rely on checkInteracion since CLL does not
                          have any notion of rigids */
                       aID=grid[ix][iy][iz].fixedID;
                       while(aID!=MAX_UINT)
                         {
                           a2=&(fixedAtoms[aID]);
                               
                           MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                           d=dx*dx+dy*dy+dz*dz;
                               
                           if (d<=d2)
                             c->push_back(AtomPair(a1,a2,d));

                           aID=fixedLL[aID];
                         }

                       /* check interaction with mobile atoms in cell */
                       aID=grid[ix][iy][iz].mobileID;
                       while(aID!=MAX_UINT)
                         {
                           a2=&(mobileAtoms[aID]);
                               
                           if (a1->id<a2->id)  /* DIFFERENCE: Avoid detecting the same interaction twice !! */
                             {
                               MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                               d=dx*dx+dy*dy+dz*dz;
                                          
                               if (d<=d2)
                                 c->push_back(AtomPair(a1,a2,d));
                             }
                           aID=mobileLL[aID];
                         }
                     }
                 }
             }
         }
    }
}


void RigidCLL::BruteForceInteractions(RigidTransform *t,list<AtomPair> *c)
{
  unsigned int i,j;
  AtomInfo *a1,*a2;
  double dx,dy,dz,d;
  
  #if (RIGIDCLL_DEBUG)
    if (cutoff==0)
      ERROR("Using BruteForceInteractions in a RigidCLL initialized to detect steric clashes");
  #endif

  if ((rID>1)||((rID==1)&&(nFixed>0)))
    {
      /* ensure the interaction table is defined */
      if (checkInteraction==NULL)
        InitInteractions();

      /* Move the atoms */
      Move(t);

      /* perform the detection */
      for(i=0;i<nMobile;i++)
        {
          a1=&(mobileAtoms[i]);

          if (checkInteraction[a1->rID][0])
            {
              for(j=0;j<nFixed;j++)
                {
                  a2=&(fixedAtoms[j]);
                     
                  MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                  d=dx*dx+dy*dy+dz*dz;
                      
                  if (d<=d2)
                    c->push_back(AtomPair(a1,a2,d));
                }
            }
        
          for(j=i+1;j<nMobile;j++)
            {
              a2=&(mobileAtoms[j]);
                     
              if (checkInteraction[a1->rID][a2->rID])
                {            
                  MANHATTAN_DISTANCE(dx,dy,dz,a1,a2);
                  d=dx*dx+dy*dy+dz*dz;
                      
                  if (d<=d2)
                    c->push_back(AtomPair(a1,a2,d));
                }
            }
        }
    }
}