Logo Search packages:      
Sourcecode: adun.app version File versions  Download package

CellListHandler.m

/*
   Project: Adun

   Copyright (C) 2005 Michael Johnston & Jordi Villa-Freixa

   Author: Michael Johnston

   Created: 2005-06-23 11:06:55 +0200 by michael johnston

   This application is free software; you can redistribute it and/or
   modify it under the terms of the GNU General Public
   License as published by the Free Software Foundation; either
   version 2 of the License, or (at your option) any later version.

   This application is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   Library General Public License for more details.

   You should have received a copy of the GNU General Public
   License along with this library; if not, write to the Free
   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111 USA.
*/
#include "Components/CellListHandler.h"

@implementation CellListHandler

/**********************************

Cell Assignment 

***********************************/

//allocates the required coordinate matrices

- (void) _initialisationForCoordinates
{
      atomCells = [memoryManager allocateIntMatrixWithRows: coordinates->no_rows withColumns: 3];
      cellNumber = [memoryManager allocateArrayOfSize: coordinates->no_rows*sizeof(int)];
}

//allocates the required interaction arrays

- (void) _initialisationForInteractions
{
      int i;

      updateCheckInteractions = [memoryManager allocateArrayOfSize: [interactions count]*sizeof(IntArrayStruct)];
      for(i = 0; i<(int)[interactions count]; i++)
            updateCheckInteractions[i].array = NULL;
}

//frees the coordinate related matrices

- (void) _clearCoordinateMatrices
{
      [memoryManager freeIntMatrix: atomCells];
      [memoryManager freeArray: cellNumber];
}     
      
//Assign each atom an x,y and z cell index 
//which then gives you the cell its in.

- (void) _initContentsArrays
{
      int i;

      if(cellContentsMatrix == NULL)
      {     
            NSDebugLLog(@"CellListHandler", @"Allocating cellContentsMatrix (%@)", NSStringFromClass([self class]));
            cellContentsMatrix = [memoryManager allocateArrayOfSize: numberOfCells*sizeof(IntArrayStruct)];
      }     
      else
            for(i=0; i< numberOfCells; i++)
                  free(cellContentsMatrix[i].array);
                        
      baseSize = (int)ceil(((double)coordinates->no_rows)/(double)numberOfCells);
      for(i=0; i<numberOfCells; i++)
      {     
            cellContentsMatrix.array = (int*)malloc(baseSize*sizeof(int));
            cellContentsMatrix.length = 0;
      }
}

//For each atom assign
//it an x,y and z cell index which then gives you the cell its in.
//For more information on the process see _updateCellIndexes:

- (void) _assignCellIndexes 
{
      int i, j;
      IntArrayStruct* cellContentsArray;

      //this exception is generic since we dont want it to be caught with the one below

      if(coordinates == NULL)
            [NSException raise: NSGenericException 
                  format: @"No atom coordinates set. Cannot build interaction list"];

      //This step is probably inefficent. Must be a better way to update the cell contents

      [self _initContentsArrays];

      NSDebugLLog(@"CellListHandler", @"Assigning atoms to cells");

      //find the cell coordinates and cell number for each atom
      //add that atom to atomsInCell array

      for(i=0; i<coordinates->no_rows; i++)
      {
            for(j=0; j<3; j++)
                  atomCells->matrix[i][j] = (int)floor((coordinates->matrix[i][j] - minSpaceBoundry.vector[j])/cellSize);
            
            cellNumber = (cellsPerAxis*cellsPerAxis)*atomCells->matrix[i][0];
            cellNumber[i] += cellsPerAxis[2]*atomCells->matrix[i][1]; 
            cellNumber[i] += atomCells->matrix[i][2];

            cellContentsArray = &cellContentsMatrix[cellNumber[i]];
            if(cellNumber[i] > numberOfCells || cellNumber[i] < 0)
                  [NSException raise: NSInternalInconsistencyException format: @"Coordinate space has changed"];

            if(cellContentsArray->length >= baseSize)
                  cellContentsArray->array = realloc(cellContentsArray->array, (cellContentsArray->length + 1)*sizeof(int));
            
            cellContentsArray->array[cellContentsArray->length] = i;
            cellContentsArray->length++; 
      }

      //trim the cell contents arrays where necessary

      for(i=0; i<numberOfCells; i++)
            cellContentsMatrix[i].array = realloc(cellContentsMatrix[i].array, cellContentsMatrix[i].length*sizeof(int));

      NSDebugLLog(@"CellListHandler", @"Complete");
}     

- (void) _updateCellIndexes 
{
      int i,j,k;
      int number;
      int *new;
      IntArrayStruct* cellContentsArray;

      //this exception is generic since we dont want it to be caught with the one below

      if(coordinates == NULL)
            [NSException raise: NSGenericException 
                  format: @"No atom coordinates set. Cannot build interaction list"];

      for(i=0; i<coordinates->no_rows; i++)
      {
            /*
              Calculate what cell this atom is now in. 
              Specfically we are calculating how many "cell blocks" it is away 
              from the origin position in each direction.

              Although we have broken the space occupied by the atoms into
              cells it is likely that atoms will move out of this space. 
              When this happens the number of cell blocks the atom is away along
              one or more axis will either be greater than the current number of cells
              we have created in that direction or it will be a negative number.

              We check if this has happened in the following loop 
              by comparing the cell index assigned to the number of cells there are
              in that direction. If it is greater or equal we have to extend the cell
              space to encompass all the atoms again. (Remember the first cell in
              each direction has index 0 following the C convention for arrays).

              We must do the same thing if the index < 0
            */

            for(j=0; j<3; j++)
            {
                  atomCells->matrix = (int)floor((coordinates->matrix[i][j] - minSpaceBoundry.vector[j])/cellSize);
                  if(atomCells->matrix[i][j] >= cellsPerAxis[j] || atomCells->matrix[i][j] < 0)
                        [NSException raise: NSInternalInconsistencyException
                               format: @"Coordinate space has changed"];
            }

            /*
              Calculate the index of the cell this atom is in.
              The cells are assigned indexes (single numbers) by counting first
              along the z axis followed by the y axis and finally the x axis 
              e.g. with two cell per axis the assignment order is
              0,0,0 -> 0,0,1 -> 0,1,0 -> 0,1,1 -> 1,0,0 -> 1,0,1, -> 1,1,0 -> 1,1,1
              The first cell will be given the index '0', the last '7'
            */

            number = (cellsPerAxis*cellsPerAxis)*atomCells->matrix[i][0];
            number += cellsPerAxis[2]*atomCells->matrix[i][1]; 
            number += atomCells->matrix[i][2];

            if(number != cellNumber[i])
            {
                  //remove from old array 
                  //\note - make these functions
                  
                  cellContentsArray = &cellContentsMatrix;
                  new = malloc((cellContentsArray->length - 1)*sizeof(int));
                  
                  for(k = 0, j=0; k< cellContentsArray->length; k++)
                        if(cellContentsArray->array[k] != i)
                        {
                              new = cellContentsArray->array;
                              j++;
                        }     

                  cellContentsArray->length = cellContentsArray->length -1;
                  free(cellContentsArray->array);
                  cellContentsArray->array = new;
                  
                  if(j != cellContentsArray->length)
                        NSWarnLog(@"Warning %d %d", j, cellContentsArray->length);

                  //add to new array
                  
                  cellContentsArray = &cellContentsMatrix;
                  new = malloc((cellContentsArray->length + 1)*sizeof(int));
                  for(k = 0, j=0 ; k < cellContentsArray->length; k++)
                  {
                        if(cellContentsArray->array[k] > i && j == k)
                        {
                              new[j] = i;
                              j++;
                        }

                        new = cellContentsArray->array;
                        j++;
                  }           

                  if(j == k)
                        new[j] = i;
                  
                  cellContentsArray->length++;
                  free(cellContentsArray->array);
                  cellContentsArray->array = new;

                  cellNumber = number;
            }
      }
}     

/**********************

Create List Elements

**********************/

- (void) _clearListDependants
{
      int i;
      
      if(updateCheckInteractions != NULL)
            for(i=0; i< (int)[interactions count]; i++)
                  free(updateCheckInteractions[i].array);
}     

- (BOOL) _checkInteractionBetweenAtomOne: (int) atomOne atomTwo: (int) atomTwo
{
      int l;
      ListElement* el_p;
      Vector3D seperation;

      for(l = 0; l < 3; l++)
            seperation.vector[l] = coordinates->matrix[atomOne][l] - coordinates->matrix[atomTwo][l];

      Ad3DVectorLengthSquared(&seperation);

      if(seperation.length < cutoff_sq)
      {
            el_p =      getElement(nonbondedList, @selector(getNewListElement));
            el_p->bond = atomOne;
            el_p->bond = atomTwo;
            if(parameters != NULL)
            {
                  el_p->params = parameters->table*parameters->table;
                  el_p->params = parameters->table*parameters->table;
            }     
            el_p->length = 0;
            return YES;
      }
            
      return NO;
}

/***********************************

Public Methods 

************************************/

- (void) createList
{
      int i,j;
      int atomIndex, cell, holder;
      int check, nocheck;
      int* inter;
      uint_fast8_t *bin; 
      int (*interact)(id, SEL, int*, int, NSRangePointer);
      ListElement* el_p;
      IntArrayStruct interactionBuffer, neighbourCells, nocheckBuffer;
      IntArrayStruct *checkInteractions, *cellContentsBuffer;
      NSEnumerator* interactionsEnum;
      NSMutableIndexSet *interaction;

      //if the cell matrices havent been created create them now

      if(!cellsInitialised)
            [self initialiseCells];

      NSDebugLLog(@"CellListHandler", @"Building Interaction List (%@).", NSStringFromClass([self class]));

      //if a list already exists free it and any other array that are allocated here

      if(listCreated)
      {
            [nonbondedList release];
            [self _clearListDependants];
      }
      
      //create a new list and deref getNewListElement
      
      nonbondedList = [AdLinkedList new];
      getElement = (ListElement* (*)(id, SEL))[nonbondedList methodForSelector:@selector(getNewListElement)];

      //New coordinates may have been set that are not accommadated
      //by the current cell space. We must catch this eventuality and rebuild
      //the cell space

      NS_DURING
      {
            [self _assignCellIndexes];
      }     
      NS_HANDLER
      {           
            NSDebugLLog(@"CellListHandler", @"The coordinates space has changed. Recalculating the cell space");
            [self clearCellMatrices];
            [self initialiseCells];
            [self _assignCellIndexes];
      }
      NS_ENDHANDLER
      
      interactionBuffer.array = (int*)malloc(coordinates->no_rows*sizeof(int));
      nocheckBuffer.array = (int*)malloc(coordinates->no_rows*sizeof(int));
      inter = interactionBuffer.array;
      bin = (uint_fast8_t*)calloc(coordinates->no_rows,sizeof(uint_fast8_t));

      atomIndex = 0;
      interactionsEnum = [interactions objectEnumerator];

      while((interaction = [interactionsEnum nextObject]))
      {
            if(([interaction count] == 0))
                        continue;
            
            //load this atoms interaction into a buffer

            interact = (int (*)(id, SEL, int*, int, NSRangePointer))[interaction methodForSelector: 
                              @selector(getIndexes:maxCount:inIndexRange:)];
            interactionBuffer.length = interact(interaction, 
                                          @selector(getIndexes:maxCount:inIndexRange:), 
                                          interactionBuffer.array, 
                                          coordinates->no_rows,
                                          NULL);
            
            for(i=0; i<interactionBuffer.length; i++)
                  bin[interactionBuffer.array[i]] = 1;
            
            //get the atoms in the current cell

            cellContentsBuffer = &cellContentsMatrix;
            for(nocheck=0, i = 0; i<cellContentsBuffer->length; i++)
            {
                  holder = cellContentsBuffer->array;
                  if(bin[holder])
                  {
                        nocheckBuffer.array = holder;
                        nocheck++;
                  }
            }

            //load the neighbourcell contents into a buffer
            
            neighbourCells = cellNeighbourMatrix;
            for(check=0, i=0; i<neighbourCells.length; i++)
            {
                  cell = neighbourCells.array;
                  cellContentsBuffer = &cellContentsMatrix;
                  if(cellContentsBuffer->length != 0)
                        for(j=0; j < cellContentsBuffer->length; j++)
                        {
                              holder = cellContentsBuffer->array[j];
                              if(bin[holder])
                              {
                                    inter[check] = holder;
                                    check++;
                              }
                        }
            }

            checkInteractions = &updateCheckInteractions;
            checkInteractions->array = malloc((check+nocheck)*sizeof(int));
            checkInteractions->length = 0;

            for(j=0; j<check; j++)
                  if([self _checkInteractionBetweenAtomOne: atomIndex atomTwo: inter[j]])
                  {
                        checkInteractions->array[checkInteractions->length] = inter;
                        checkInteractions->length++;
                  }                             
      

            for(i=0;i<nocheck; i++)
            {     
                  el_p =      getElement(nonbondedList, @selector(getNewListElement));
                  el_p->bond = atomIndex;
                  el_p->bond = nocheckBuffer.array;
                  if(parameters != NULL)
                  {
                        el_p->params = parameters->table*parameters->table;
                        el_p->params = parameters->table*parameters->table;
                  }
                  el_p->length = 0;
                  checkInteractions->array[checkInteractions->length] = nocheckBuffer.array;
                  checkInteractions->length++;
            }
            
            checkInteractions->array = realloc(checkInteractions->array, checkInteractions->length*sizeof(int));
            memset(bin, 0, coordinates->no_rows*sizeof(uint_fast8_t));
            atomIndex++;
      }

      GSPrintf(stderr, @"There are %d nonbonded interactions.\n",[nonbondedList listCount]);
      free(interactionBuffer.array);
      free(nocheckBuffer.array);
      free(bin);
      listCreated = YES;
}     

/*******************************

Update and related functions

********************************/

- (void) _updateListRemovingTo: (IntArrayStruct*) removedInteractions
{
      int atomOne, atomTwo;
      ListElement* holder, *list_p;

      list_p = [nonbondedList listStart]->next;
      while(list_p->next != NULL)
      {     
            if(list_p->length > cutoff)
            {
                  holder = list_p;
                  list_p = list_p->next;

                  atomOne = holder->bond;
                  atomTwo = holder->bond;
                  
                  if(removedInteractions[atomOne].length >= updateCheckInteractions[atomOne].length)
                        removedInteractions.array = realloc(removedInteractions[atomOne].array, 
                                                (removedInteractions[atomOne].length + 1)*sizeof(int));
                  
                  removedInteractions.array[removedInteractions.length] = atomTwo;  
                  removedInteractions.length++;

                  [nonbondedList freeListElement: holder];
            }     
            else
                  list_p = list_p->next;
      }     
}

/*
The update algorithm
1) Update all the cells contents
2) Remove all interactions beyond the cutoff from the list
3) Add any interactions that are now within the cutoff to the list
*/

- (void) update
{
      int i, j, k, cell, holder;
      int atomTwo;
      IntArrayStruct currentInteractions, neighbourCells, *checkInteractions, newInteractions;
      IntArrayStruct *cellContentsBuffer, *removedInteractions, *removed;
      uint_fast8_t *cellBin;
      NSMutableIndexSet* interaction;
      Vector3D separation;
      int (*interact)(id, SEL, int);      
      id (*fetch)(id, SEL, int);
      SEL selector, selector2;

      NSDebugLLog(@"CellListHandler", @"Updating Lists");

      /*
         * Setup
       */

      cellBin = (uint_fast8_t*)calloc(sizeof(uint_fast8_t),coordinates->no_rows);
      currentInteractions.array = malloc(coordinates->no_rows*sizeof(int));
      removedInteractions = [memoryManager allocateArrayOfSize: 
                        [interactions count]*sizeof(IntArrayStruct)];
      
      selector = @selector(containsIndex:);
      selector2 = @selector(objectAtIndex:);

      //we cache this method call since will be using it alot
      fetch = (id (*)(id, SEL, int))[interactions 
                  methodForSelector: @selector(objectAtIndex:)];

      for(i=0; i<(int)[interactions count]; i++)
      {
            removedInteractions.array = (int*)malloc(updateCheckInteractions[i].length*sizeof(int));
            removedInteractions.length = 0;
      }
      
      /*
         * Step 1. Update cell contents
         */

      NSDebugLLog(@"CellListHandler", @"Assigning atoms to cells");
      
      NS_DURING
      {
            [self _updateCellIndexes];
      }     
      NS_HANDLER
      {           
            NSDebugLLog(@"CellListHandler",
                   @"The coordinates space has changed. Recalculating the cell space");
            [self clearCellMatrices];
            [self initialiseCells];
            [self _assignCellIndexes];
      }
      NS_ENDHANDLER

      NSDebugLLog(@"CellListHandler", @"Building List");
      
      /*
       * Step 2. Remove interactions greater than cutoff from list
       */
      
      [self _updateListRemovingTo: removedInteractions];
      
      /*
         Step 3. Add new interactions.
         This is the most complicated step since -
         A) We dont want to check the same interactions that were
            checked in step 2 i.e. those that are still in the list
            and those we removed.
         B) We dont want to check more atoms than we have to.

         We check interactions atom by atom starting with the first.
      */

      holder = cutoff + diagonal;
      for(k=0; k< (int)[interactions count] - 1; k++)
      {     
            interaction = fetch(interactions, selector2, k);
            if([interaction count] == 0)
                  continue;

            currentInteractions.length = 0;
            interact = (int (*)(id, SEL, int))[interaction methodForSelector: selector];
            
            /*
              Load the possible interactions into a buffer.  This is all
              other atoms in this atoms cell plus all atoms in neigbouring 
              cells (aslong as part of the neighbouring cell is within the cutoff). 

              A part of a neighbouring cell is within the cutoff distance of the
              atom in question if the distance from the atom to the center of the
              neighbouring cell is less than the distance defined by the variable holder.
            */

            //get the atoms in the current cell
            cellContentsBuffer = &cellContentsMatrix;
            for(i = 0; i<cellContentsBuffer->length; i++)
                  if(cellContentsBuffer->array[i] > k)
                  {
                        currentInteractions.array[currentInteractions.length] = cellContentsBuffer->array;
                        currentInteractions.length++;
                  }

            //load the neighbourcell contents 
            neighbourCells = cellNeighbourMatrix;
            for(i=0; i<neighbourCells.length; i++)
            {
                  cell = neighbourCells.array;
                  cellContentsBuffer = &cellContentsMatrix;
                  if(cellContentsBuffer->length != 0 && 
                        cellContentsBuffer->array[cellContentsBuffer->length-1] > k)
                  {
                        for(j=0; j<3; j++)
                              separation.vector[j] = cellCenterMatrix->matrix[cell][j] -
                                                 coordinates->matrix[k][j];
                        
                        Ad3DVectorLength(&separation);
                        //If the separation is less than holder than some of this
                        //cells atoms could be within the cutoff
                        if(separation.length < holder) 
                              for(j=0; j < cellContentsBuffer->length; j++)
                                    if(cellContentsBuffer->array[j] > k)
                                    {
                                          currentInteractions.array[currentInteractions.length] = 
                                                      cellContentsBuffer->array;
                                          currentInteractions.length++;
                                    }
                  }
            }
            
            /*
              As mentioned above we dont want to check any interactions that were 
              in the list when we began updating since these were delt with in step 2.
              Each entry in updateCheckInteractions is an array that contains the indexes 
              of all the atoms that this atom was interacting with in the last list i.e.
              the ones we dont want to check. 

              In order to avoid these use a "bin". This is an array with one element
              for every atom. The element is set to 1 if the interaction was in the last
              list and is zero otherwise.
            */

            //Set the relevant bin elements to 1
            checkInteractions = &updateCheckInteractions;
            for(i=0; i<checkInteractions->length; i++)
                  cellBin[checkInteractions->array[i]] = 1;

            /*
              Go through all the possible interactions we added to the buffer above
              (currentInteractions). Skip those whose element has been set to one 
              in the bin and any who arent in this atoms interaction list.
            */

            newInteractions.length = 0;
            newInteractions.array = malloc(coordinates->no_rows*sizeof(int));
            for(i=0; i< currentInteractions.length; i++)
            {
                  atomTwo = currentInteractions.array;
                  if(cellBin[atomTwo] != 1)
                        if(interact(interaction, selector, atomTwo))
                              if([self _checkInteractionBetweenAtomOne: k atomTwo: atomTwo])
                              {
                                    newInteractions.array[newInteractions.length] = atomTwo;
                                    newInteractions.length++;
                              }
            }

            /*
              Finally we want to create a new updateCheckInteraction entry for this
              atom. Therefore we need to remove all the atoms that we removed in 
              step 2. 

              removedInteractions contains the indexes of the atoms we removed. We
              set the corresponding element in the bin to 0.

              After this the only elements in the bin set to one are those
              that remained in the list after step 2. These combined with the ones
              we added above (in newInteractions) are what we want.
            */

            removed = &removedInteractions;
            for(i=0; i<removed->length; i++)
                  cellBin[removed->array[i]] = 0;

            //add the interactions still in the list to newInteractions

            for(i=0; i<checkInteractions->length; i++)      
                  if(cellBin[checkInteractions->array[i]] == 1)
                  {
                        newInteractions.array[newInteractions.length] = checkInteractions->array;
                        newInteractions.length++;
                  }                       

            newInteractions.array = realloc(newInteractions.array, newInteractions.length*sizeof(int));
            free(checkInteractions->array);
            free(removed->array);
            checkInteractions->array = newInteractions.array;
            checkInteractions->length = newInteractions.length;
            memset(cellBin, 0, sizeof(uint_fast8_t)*coordinates->no_rows);
      }

      NSDebugLLog(@"CellListHandler", @"Update Complete");

      free(cellBin);
      free(currentInteractions.array);
      free(removedInteractions);
}

/*********************************

Object Creation and Maintenance Methods 

**********************************/

- (void) _useEnvironmentDefaults
{
      cutoff = 12;
      maxSpaceSize = 100;
}

- (void) _initialiseDependants
{
      cutoff_sq = cutoff*cutoff;
      cellSize = cutoff/2;
      diagonal = 0.5*sqrt(2)*cellSize;
      cellCut = diagonal + cutoff;
      inCut = cutoff/sqrt(2);
      cellCut *= cellCut;
}

- (id) init
{
      return [self initWithEnvironment: nil];
}

00720 - (id) initWithEnvironment: (id) object
{
      return [self initWithEnvironment: object observe: YES];
}

00725 - (id) initWithEnvironment: (id) object observe: (BOOL) value
{
      if((self = [super initWithEnvironment: object observe: value]))
      {
            parameters = NULL;
            memoryManager = [AdMemoryManager appMemoryManager];
            updateCheckInteractions = NULL;
            if(environment == nil)
                  [self _useEnvironmentDefaults];
            else
            {
                  [self synchroniseWithEnvironment];
                  [self registerWithEnvironment];
            }

            [self _initialiseDependants];
            
      }

      return self;
}

- (void) dealloc
{
      [self clearCellMatrices];
      [self _clearCoordinateMatrices];
      [nonbondedList release];
      [self _clearListDependants];
      free(cellsPerAxis);
      free(updateCheckInteractions);
      [super dealloc];
}

/*
 * Environment observation
 */

00762 - (void) updateForKey: (NSString*) key value: (id) value object: (id) object
{
      if([key isEqual: @"Cutoff"])
      {
            cutoff = [value floatValue];
            [self _initialiseDependants]; 
      }
}

00771 - (void) registerWithEnvironment
{
      if(observesEnvironment)
            [environment addObserver: self forKey: @"CutOff"];
}

00777 - (void) deregisterWithEnvironment
{
      [environment removeObserver: self forKey: @"CutOff"];
}

00782 - (void) synchroniseWithEnvironment
{
      cutoff = [[environment valueForKey: @"CutOff"] 
                  floatValue];

      maxSpaceSize = [[environment valueForKey: @"MaximumSpaceSize"]
                        doubleValue]; 

      //check maxSpaceSize exists - we should really do this for all 
      //environment variables e.g. Record names of variable the
      //enviroment didnt hold and use defaults

      if(maxSpaceSize == 0)
            maxSpaceSize = 100;
}

00798 - (void) setEnvironment: (id) object
{
      [self deregisterWithEnvironment];
      object = environment;
      [self registerWithEnvironment];
}


/*******************

Accessors

*********************/

- (void) setCoordinates: (NSValue*) aValue
{
      if(coordinates != NULL)
            [self _clearCoordinateMatrices];
                  
      coordinates = [aValue pointerValue];
      [self _initialisationForCoordinates];
}

- (void) setNonbondedTopology: (NSArray*) anArray
{
      [self _clearListDependants];
      if(updateCheckInteractions != NULL)
            free(updateCheckInteractions);

      interactions = anArray; 
      [self _initialisationForInteractions];
}

- (void) setParameters: (NSValue*) aValue
{
      parameters = [aValue pointerValue];
}

- (void) setCutoff: (int) aValue
{
      cutoff = aValue;
      [self _initialiseDependants];
}

- (NSValue*) nonbondedInteractions
{
      return [NSValue valueWithPointer: [nonbondedList listStart]];
}

- (int) numberNonbondedInteractions
{
      return [nonbondedList listCount];
}

/**********************

Coding

**********************/

- (id) initWithCoder: (NSCoder*) decoder
{
      self = [super initWithCoder: decoder];
      if([decoder allowsKeyedCoding])
      {
            memoryManager = [AdMemoryManager appMemoryManager];
            environment = [AdEnvironment globalEnvironment];
            if(environment != nil)
            {
                  [self synchroniseWithEnvironment];
                  [self registerWithEnvironment];
            }
            else
                  [self _useEnvironmentDefaults];
            
            [self _initialiseDependants];
            parameters = NULL;
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"%@ class does not support non keyed coding", [self class]];

      return self;
}

- (void) encodeWithCoder: (NSCoder*) encoder
{
      [super encodeWithCoder: encoder];
}

@end


/**
Preliminary refactoring of cell related methods
to a category. The Cell assignment methods are
still in the main class but need at bit of changing
before they can be moved here
*/

@implementation CellListHandler (CellMaintainence)

/**
Locates the boundaries of the rectangular box that encloses all the atoms. 
We then extend them so we can fit an integer number of cells in each 
direction.  (The cell size is half the cutoff chosen).
Finally we add some padding by extending each boundary by half the cell size 
- Effectivly adding one extra cell in each direction.

To handle exploding simulations we check against the user defined max
space size. If the cell space would be greater in any direction than
this number we raise an exception.
*/

00912 - (void) _locateCellSpaceBoundries
{
      int i;
      double* coordinateArray;
      double xdiff, ydiff, zdiff;

      coordinateArray = [memoryManager allocateArrayOfSize: coordinates->no_rows*sizeof(double)];

      NSDebugLLog(@"CellListHandler", @"Finding the space extremes");
      
      //the mininum and maximum x coordinates   

      for(i=0; i<coordinates->no_rows; i++)
            coordinateArray[i] = coordinates->matrix[i][0];

      minSpaceBoundry.vector[0] = coordinateArray[AdDoubleArrayMin(coordinateArray, coordinates->no_rows)];
      maxSpaceBoundry.vector = coordinateArray;

      //the mininum and maximum y coordinates   

      for(i=0; i<coordinates->no_rows; i++)
            coordinateArray[i] = coordinates->matrix[i][1];

      minSpaceBoundry.vector[1] = coordinateArray[AdDoubleArrayMin(coordinateArray, coordinates->no_rows)];
      maxSpaceBoundry.vector = coordinateArray;
      
      //the mininum and maximum z coordinates   

      for(i=0; i<coordinates->no_rows; i++)
            coordinateArray[i] = coordinates->matrix[i][2];

      minSpaceBoundry.vector[2] = coordinateArray[AdDoubleArrayMin(coordinateArray, coordinates->no_rows)];
      maxSpaceBoundry.vector = coordinateArray;

      free(coordinateArray);

      NSDebugLLog(@"CellListHandler", @"Minimum: %-10.3lf %-10.3lf %-10.3lf", 
                  minSpaceBoundry.vector[0], minSpaceBoundry.vector[1], minSpaceBoundry.vector[2]);
      NSDebugLLog(@"CellListHandler", @"Maximim: %-10.3lf %-10.3lf %-10.3lf", 
                  maxSpaceBoundry.vector[0], maxSpaceBoundry.vector[1], maxSpaceBoundry.vector[2]);

      //find how many cells are needed in each direction
      
      cellSpaceDimensions.vector = maxSpaceBoundry.vector - minSpaceBoundry.vector;
      cellSpaceDimensions.vector = maxSpaceBoundry.vector - minSpaceBoundry.vector;
      cellSpaceDimensions.vector = maxSpaceBoundry.vector - minSpaceBoundry.vector;

      //if by chance any of the cellSpaceDimensions == 0, put it equal to one
      //this will happen if the molecule lies entirely in an axis plane.

      for(i=0; i<3; i++)
            if(cellSpaceDimensions.vector[i] == 0)
                  cellSpaceDimensions.vector = 1;

      cellsPerAxis = [memoryManager allocateArrayOfSize: 3*sizeof(int)];
      cellsPerAxis = (int)ceil(cellSpaceDimensions.vector[0]/cellSize);
      cellsPerAxis = (int)ceil(cellSpaceDimensions.vector[1]/cellSize);
      cellsPerAxis = (int)ceil(cellSpaceDimensions.vector[2]/cellSize);

      NSDebugLLog(@"CellListHandler", @"Cells per axis: %d %d %d", 
                  cellsPerAxis[0], 
                  cellsPerAxis[1],
                  cellsPerAxis[2]);

      //calculate how much we need to move each boundary to accomadate the necessary cells

      xdiff = cellsPerAxis*cellSize - cellSpaceDimensions.vector;
      ydiff = cellsPerAxis*cellSize - cellSpaceDimensions.vector;
      zdiff = cellsPerAxis*cellSize - cellSpaceDimensions.vector;

      //increase each boundary by half the difference

      minSpaceBoundry.vector -= xdiff/2;
      minSpaceBoundry.vector -= ydiff/2;
      minSpaceBoundry.vector -= zdiff/2;
      maxSpaceBoundry.vector += xdiff/2;
      maxSpaceBoundry.vector += ydiff/2;
      maxSpaceBoundry.vector += zdiff/2;
      
      NSDebugLLog(@"CellListHandler", @"Recalculated  space extremes");
      NSDebugLLog(@"CellListHandler", @"Minimum: %-10.3lf %-10.3lf %-10.3lf", 
                  minSpaceBoundry.vector[0], minSpaceBoundry.vector[1], minSpaceBoundry.vector[2]);
      NSDebugLLog(@"CellListHandler", @"Maximim: %-10.3lf %-10.3lf %-10.3lf", 
                  maxSpaceBoundry.vector[0], maxSpaceBoundry.vector[1], maxSpaceBoundry.vector[2]);

      //recalculate the cellSpaceDimensions
      //check that none of the cellSpaceDimensions are greater the maxSpaceSize

      for(i=0; i<3; i++)
      {
            cellSpaceDimensions.vector = cellsPerAxis*cellSize;
            if(cellSpaceDimensions.vector[i] > maxSpaceSize)
                  [NSException raise: NSInternalInconsistencyException
                         format: @"Coordinate space has exceeded size restriction (%lf, %lf)",
                         maxSpaceSize, cellSpaceDimensions.vector];
      }

      NSDebugLLog(@"CellListHandler", @"The cell space dimensions are %lf, %lf, %lf", 
                  cellSpaceDimensions.vector[0],
                  cellSpaceDimensions.vector[1], 
                  cellSpaceDimensions.vector[2]);
}

//we first need to divide the space into a series of cells
//with defined centers and corresponding indexes (0,0,0) (0,0,1) etc.

- (void) _createCellMatrices
{
      int i,j,k;
      int xIndex, yIndex, zIndex, currentCell, number;
      IntArrayStruct* cellNeighbours;

      for(numberOfCells=1, i=0; i<3; i++)
            numberOfCells *= cellsPerAxis[i]; 

      NSDebugLLog(@"CellListHandler", @"Creating the cell matrices. There are %d cells", numberOfCells);
      
      cellCenterMatrix = [memoryManager allocateMatrixWithRows: numberOfCells withColumns: 3];
      cellIndexMatrix = [memoryManager allocateMatrixWithRows: numberOfCells withColumns: 3];

      //find center of cell (0,0,0) (the origin cell)

      originCellCenter.vector = minSpaceBoundry.vector + cellSize/2;
      originCellCenter.vector = minSpaceBoundry.vector + cellSize/2;
      originCellCenter.vector = minSpaceBoundry.vector + cellSize/2;
      
      NSDebugLLog(@"CellListHandler", @"The origin cell center is %lf, %lf, %lf", 
                  originCellCenter.vector[0], originCellCenter.vector[1], originCellCenter.vector[2]);

      //the cell order is (x,y,z) - assign cell indexes in this order

      for(currentCell=0, i=0; i<cellsPerAxis;i++)
            for(j=0; j<cellsPerAxis;j++)
                  for(k=0; k<cellsPerAxis;k++)
                  {
                        cellIndexMatrix->matrix = i;
                        cellIndexMatrix->matrix = j;
                        cellIndexMatrix->matrix = k;
                        currentCell++;
                  }
      
      //use the cellIndexMatrix to assign coordinates to the cell centers

      for(i=0; i<numberOfCells; i++)
            for(j=0; j<3; j++)
                  cellCenterMatrix->matrix[i][j] = originCellCenter.vector[j] + cellIndexMatrix->matrix[i][j]*cellSize;

      //create the cell neighbour array

      cellNeighbourMatrix = [memoryManager allocateArrayOfSize: numberOfCells*sizeof(IntArrayStruct)];
      for(currentCell = 0; currentCell < numberOfCells; currentCell++)
      {
            cellNeighbours = &cellNeighbourMatrix;
            cellNeighbours->array = (int*)malloc(80*sizeof(int));
            cellNeighbours->length = 0;

            xIndex = (int)cellIndexMatrix->matrix[currentCell][0];
            yIndex = (int)cellIndexMatrix->matrix[currentCell][1];
            zIndex = (int)cellIndexMatrix->matrix[currentCell][2];

            for(i= xIndex - 2; i <= xIndex + 2; i++)
                  if(i >= 0 && i < cellsPerAxis[0])
                        for(j = yIndex -2; j<= yIndex + 2; j++)
                              if(j >= 0 && j < cellsPerAxis[1])
                                    for(k = zIndex -2; k <= zIndex + 2; k++)
                                          if(k>=0 && k < cellsPerAxis[2])
                                          {
                                                number =  (cellsPerAxis[2]*cellsPerAxis[1])*i + cellsPerAxis[2]*j + k;
                                                if(number != currentCell && number < numberOfCells)
                                                {
                                                      //check if reallocing is needed
                                                      if(cellNeighbourMatrix[currentCell].length >= 80)
                                                            cellNeighbours->array = realloc(cellNeighbours->array,
                                                                         (cellNeighbours->length + 1)*sizeof(int));
                                                      cellNeighbours->array[cellNeighbours->length] = number;
                                                      cellNeighbours->length++;
            
                                                }
                                          }

            //trim the memory to the neccessary size

            cellNeighbours->array = realloc(cellNeighbours->array, (cellNeighbours->length)*sizeof(int));
      }     

      //initialise CellContents array

      cellContentsMatrix = [memoryManager allocateArrayOfSize: numberOfCells*sizeof(IntArrayStruct)];
      baseSize = (int)ceil(((double)coordinates->no_rows)/(double)numberOfCells);

      for(i=0; i<numberOfCells; i++)
      {     
            cellContentsMatrix.array = (int*)malloc(baseSize*sizeof(int));
            cellContentsMatrix.length = 0;
      }
}

- (void) initialiseCells
{
      NSDebugLLog(@"CellListHandler", @"Initialising cell space");
      
      [self _locateCellSpaceBoundries];
      [self _createCellMatrices];

      cellsInitialised = YES;
}

- (void) clearCellMatrices
{
      int i;

      [memoryManager freeMatrix: cellCenterMatrix];
      [memoryManager freeMatrix: cellIndexMatrix];
      
      for(i=0; i<numberOfCells;i++)
      {
            free(cellNeighbourMatrix[i].array);
            free(cellContentsMatrix[i].array);
      }
      free(cellNeighbourMatrix);
      free(cellContentsMatrix);
}

@end

Generated by  Doxygen 1.6.0   Back to index