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

SphericalBox.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/SphericalBox.h"

@implementation SphericalBox

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

Expanding the topology 

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

-(InterTable*) _expandTable:(InterTable*) table 
            times: (const int) expansion indexColumns:(const int) number
{
      int i, j, k, index, exp_rows, solute_offset;
      InterTable *expanded_table;
      
      //malloc memory for the expanded matrix
      exp_rows = table->no_interactions*expansion;
      expanded_table = [memoryManager allocateInterTableWithRows: exp_rows
                              withColumns: table->no_columns]; 

      //copy the matrix "expansion" times
      
      for(i=0; i < exp_rows; i = i + table->no_interactions)
            for(index = 0, j=i; j < i + table->no_interactions; j++, index++)
                  for(k=0; k < table->no_columns; k++)
                        expanded_table->table[j][k] = table->table[index][k];
      
      //give the correct indices to each entry - currently we just have the same indices repeated expansion times
      
      if(number!=0)
            for(i=0, index = 0; i < exp_rows; i += table->no_interactions, index += atomsPerMolecule)
                  for(j=i; j < i + table->no_interactions; j++)
                        for(k=0; k < number; k++)                 
                              expanded_table->table[j][k] += index;
      
      return expanded_table;
}     

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

Data Retrieval

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

- (void) _retrieveCoordinates
{
      int i, k, j,currentAtom;
      AdMatrix* dataSourceCoordinates;

      dataSourceCoordinates = [[dataSource objectValueForCoordinates: self] pointerValue];

      //find the solvent subsection mass and hence the number of molecules in the sphere
      
      for(solventMass=0, i=0; i<dataSourceCoordinates->no_rows; i++)
            solventMass += dataSourceCoordinates->matrix[i][3]; 

      NSDebugLLog(@"SphericalBox", @"Mass %lf", solventMass); 

      totalNumberOfMolecules = solventDensity*sphereVolume/solventMass;
      totalNumberOfAtoms = totalNumberOfMolecules*dataSourceCoordinates->no_rows;
      atomsPerMolecule = dataSourceCoordinates->no_rows;
      no_obscured_molecules = 0;
      
      NSDebugLLog(@"SphericalBox", @"Retrieving coordinates and expanding. There are %d molecules in the sphere", 
                  totalNumberOfMolecules);
      NSDebugLLog(@"SphericalBox", @"There are %d atomsPerMolecule", atomsPerMolecule);

      //expand the dataSourceAtomTypes and create a coordinate matrix

      solventCoordinates = [memoryManager allocateMatrixWithRows: totalNumberOfAtoms withColumns: 6];

      for(currentAtom = 0, i=0; i<totalNumberOfMolecules; i++)
            for(j=0;j<atomsPerMolecule; j++)
            {
                  for(k=0; k<6; k++)
                        solventCoordinates->matrix[currentAtom][k] = dataSourceCoordinates->matrix[j][k];

                  currentAtom++;
            }

      NSDebugLLog(@"SphericalBox", @"Retrieved Coordinates and Atom types");

}

//FIXME: this is hacky since we have to know the number of atoms involved in each interaction type
//but this information isnt available. 

- (void) _retrieveBondedInteractions
{
      int indexColumns = 0;
      InterTable* interactionTable;
      id interactionType, interaction, dataSourceBondedInteractions;
      NSEnumerator* interactionEnum;
      InterTable* temp;

      NSDebugLLog(@"SphericalBox", @"Retrieveing bonded interactions from data source");  

      solventBondedInteractions = [[NSMutableDictionary dictionaryWithCapacity: 1] retain];
      dataSourceBondedInteractions = [dataSource objectValueForBondedInteractions: self];
      interactionEnum = [dataSourceBondedInteractions keyEnumerator];
      while(interactionType = [interactionEnum nextObject])
      {
            if([interactionType isEqual: @"HarmonicBond"])
                  indexColumns = 2;
            else if([interactionType isEqual: @"HarmonicAngle"])
                  indexColumns = 3;
            else if([interactionType isEqual: @"FourierTorsion"])
                  indexColumns = 4;
            else if([interactionType isEqual: @"HarmonicImproperTorsion"])
                  indexColumns = 4;
            else
                  [NSException raise: NSInvalidArgumentException 
                        format: @"Unknown interaction type %@", interactionType];

            interaction = [dataSourceBondedInteractions valueForKey: interactionType];
            temp = [interaction pointerValue];
            interactionTable = [self _expandTable: temp 
                              times: currentNumberOfMolecules
                              indexColumns: indexColumns];
            [solventBondedInteractions setObject: [NSValue valueWithPointer: interactionTable]
                  forKey: interactionType];
      }
}

- (void) _retrieveNonbondedInteractions
{
      int i, j, currentAtom, index, offsetIndex;
      id dataSourceNonbonded, dataSourceNonbondedTypes, type, interaction;
      NSMutableIndexSet* indexes, *sourceIndexes, *currentInteractions;
      NSRange indexRange;
      NSEnumerator* typeEnum;
      InterTable* interactionTable;

      NSDebugLLog(@"SphericalBox", @"Retrieving Nonbonded interactions from data source");

      dataSourceNonbonded = [dataSource objectValueForNonbondedInteractions: self];
      dataSourceNonbondedTypes = [dataSource objectValueForNonbondedInteractionTypes: self];
      solventNonbonded = [[NSMutableArray arrayWithCapacity: (totalNumberOfAtoms - 1)] retain];
      solventNonbondedTypes = [[NSMutableDictionary dictionaryWithCapacity: 1] retain];

      //first create the inter molecule nonbonded list
      //we create entries for the last molecule so we can easily add an intra-molecule
      //interactions in the next step

      for(i=0; i<currentNumberOfMolecules; i++)
            for(j=0; j<atomsPerMolecule; j++)
            {
                  indexRange.location = (i+1)*atomsPerMolecule;
                  indexRange.length = currentNumberOfAtoms - indexRange.location;
                  indexes = [NSMutableIndexSet indexSetWithIndexesInRange: indexRange];
                  [solventNonbonded addObject: indexes];
            }

      //add the intra molecule nonbonded interactions
      for(currentAtom = 0,i=0; i<currentNumberOfMolecules; i++)
      {
            for(j=0; j<atomsPerMolecule - 1; j++)
            {
                  currentInteractions = [solventNonbonded objectAtIndex: currentAtom];
                  sourceIndexes = [dataSourceNonbonded objectAtIndex: j];
                  index = [sourceIndexes firstIndex];
                  while(index != NSNotFound)
                  {
                        offsetIndex = index + i*atomsPerMolecule;
                        [currentInteractions addIndex: offsetIndex];
                        index = [sourceIndexes indexGreaterThanIndex: index];
                  }

                  currentAtom++;
            }

            //skip the last atom in each molecule
            currentAtom++;
      }

      //remove the last atom entry since its not needed

      [solventNonbonded removeLastObject];

      //expand the vdw parameters

      typeEnum = [dataSourceNonbondedTypes keyEnumerator];
      while(type = [typeEnum nextObject])
      {
            interaction = [dataSourceNonbondedTypes valueForKey:type];
            if(![interaction isKindOfClass: [NSNull class]])
            {
                  interactionTable = [self _expandTable: [interaction pointerValue] 
                                    times: currentNumberOfMolecules
                                    indexColumns: 0];
                  [solventNonbondedTypes setObject: [NSValue valueWithPointer: interactionTable]
                        forKey: type];
            }
            else
                  [solventNonbondedTypes setObject: [NSNull null]
                        forKey: type];
                  
      }

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

- (void) _clearCurrentSystem
{
      NSEnumerator* anEnum;
      id obj;
      InterTable* temp;

      [solventAtomTypes release];   
      [memoryManager freeMatrix: currentSolventCoordinates];

      anEnum = [solventBondedInteractions objectEnumerator];      
      while(obj = [anEnum nextObject])
      {
            //hack for topology types with no interactions
            temp = [obj pointerValue];
            if(temp->no_interactions != 0)
                  [memoryManager freeInterTable: temp];
            else
                  free(temp);
      }
      [solventBondedInteractions release];      //clear the interactions
      
      [solventNonbonded release];
      anEnum = [solventNonbondedTypes objectEnumerator];    

      while(obj = [anEnum nextObject])
            if(![obj isKindOfClass: [NSNull class]])
            {
                  temp = [obj pointerValue];
                  [memoryManager freeInterTable: temp];
            }     

      [solventNonbondedTypes release];

}

- (void) _setCurrentSystem
{
      int i, j;
      id dataSourceAtomTypes;
      
      //use the information in currentSolventCoordinates to build system

      dataSourceAtomTypes = [dataSource objectValueForAtomTypes: self];
      solventAtomTypes = [[NSMutableArray arrayWithCapacity: currentNumberOfAtoms] retain];
      for(i=0; i<currentNumberOfMolecules; i++)
            [solventAtomTypes addObjectsFromArray: dataSourceAtomTypes];
      
      [self _retrieveBondedInteractions];
      [self _retrieveNonbondedInteractions];

      GSPrintf(stderr, @"There are %d current solvent atoms\n", currentSolventCoordinates->no_rows);
}

/**
When reloadData is called we need to free any objects associated with
a previous data load
*/

00285 - (void) _cleanUp
{
      //if solventAtomTypes is not nil then data was previously loaded

      if(solventAtomTypes != nil)
      {
            //FIXME: This code frees currentSolventCoordinates twice
            //(Once in _clearCurrentSystem). However the core has
            //never crashed from it! Trying to find out why.

            [self _clearCurrentSystem];
            [memoryManager freeMatrix: currentSolventCoordinates];
            [memoryManager freeMatrix: solventCoordinates];
      }
}

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

Populating the Box 

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

- (void) _randomlyOrientMolecules
{
      int i, k,j, currentAtom;
      double angle;
      double rotated;
      
      for(currentAtom=0, j=0; j< totalNumberOfMolecules; j++)
      {
            //generate 3 random angles between -pi and pi
            for(i=0; i<3; i++)
                  angle[i] = (2*gsl_rng_uniform(twister) -1)*M_PI;                  
      
            //rotate each molecule by the angles in angle
            //this involves rotating each atom of the molecule seperately

            for(i=0; i < atomsPerMolecule; i++)
            {     
                  AdRotate3DVector(solventCoordinates->matrix[currentAtom], angle, rotated);
                  for(k=0; k<3; k++)
                        solventCoordinates->matrix[currentAtom][k] = rotated[k];

                  currentAtom++;
            }
      }
}

- (void) _placeMoleculesAtSites: (int*) sites
{
      int i, j, k;
      int current_site, current_atom;
      AdMatrix* gridMatrix;

      //randomly orientate the molecules

      [self _randomlyOrientMolecules];

      gridMatrix = [[solventGrid valueForKey: @"grid"] pointerValue];
      for(i=0; i<totalNumberOfMolecules; i++)
      {
            current_site = sites;
            current_atom = i*atomsPerMolecule;
            for(j=current_atom; j<current_atom + atomsPerMolecule; j++)
                  for(k=0; k<3; k++)
                        solventCoordinates->matrix[j][k] += gridMatrix->matrix[current_site][k];
      }
}

- (int*) _chooseGridPoints
{
      int i, j;
      int *array, *point_array;
      AdMatrix* gridMatrix;   

      NSDebugLLog(@"SphericalBox", @"Choosing grid_points for %d molecules.\n", totalNumberOfMolecules);

      //create an array of size grid_points, containing the numbers between 0 and grid_points
      
      gridMatrix = [[solventGrid valueForKey:@"grid"] pointerValue];
      array = (int*)[memoryManager allocateArrayOfSize: gridMatrix->no_rows*sizeof(int)];
      for(i=0; i < gridMatrix->no_rows; i++)
            array[i] = i;
      
      gsl_ran_shuffle(twister, array, gridMatrix->no_rows, sizeof(int));
      
      //malloc memory for the array to hold the choosen grid points
      
      point_array = (int*)[memoryManager allocateArrayOfSize: totalNumberOfMolecules*sizeof(int)];
      for(i=0; i<totalNumberOfMolecules; i++)
            point_array[i] = array[i];

      [memoryManager freeArray: array];
      
      return point_array;     
}     

- (void) _populateBox
{
      int* points;      
      int i, j;

      //choose a number of grid points equal to the number of 
      //molecles in the box

      points = [self _chooseGridPoints];

      //place solvent molecules at each point

      [self _placeMoleculesAtSites: points];    

      //we keep a seperate matrix for coordinates containing only unobscured molecules and we use
      //this to build all the bonded and nonbonded interactions etc. The main matrix to enable easy
      //expansion and contraction of exlusion areas.

      currentSolventCoordinates = [memoryManager allocateMatrixWithRows: totalNumberOfAtoms withColumns: 6];
      for(i=0; i<totalNumberOfAtoms; i++)
            for(j=0; j<6; j++)
                  currentSolventCoordinates->matrix[i][j] = solventCoordinates->matrix[i][j];
      
      currentNumberOfAtoms = currentSolventCoordinates->no_rows;
      currentNumberOfMolecules = currentNumberOfAtoms/atomsPerMolecule;
      [memoryManager freeArray: points];
}

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

Creation and Maintainence

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

- (void) _createGrid
{
      NSMutableArray* spacing;
      
      //Create the grid on which we'll place the solvent    

      spacing = [NSMutableArray arrayWithCapacity: 1];

      //FIXME: calculate this correctly!!!!

      [spacing addObject: [NSNumber numberWithDouble: 3.0]];
      [spacing addObject: [NSNumber numberWithDouble: 3.0]];
      [spacing addObject: [NSNumber numberWithDouble: 3.0]];
      
      NSDebugLLog(@"SphericalBox", @"Grid spacing %@", spacing);

      solventGrid = [AdGrid gridWithSpacing: spacing 
                        cavity: self
                        environment: environment];
      [solventGrid retain];
}

- (void) _initialiseDependants
{
      NSMutableArray* spacing;
      id obj;

      sphereVolume = 4*M_PI*pow(sphereRadius, 3)/3;
      
      obj = [NSArray arrayWithObjects:
                  [NSNumber numberWithDouble: sphereRadius],      
                  [NSNumber numberWithDouble: -sphereRadius],
                  nil];       

      sphereExtremes = [NSArray arrayWithObjects: obj, [obj copy], [obj copy], nil];
      [sphereExtremes retain];

      NSDebugLLog(@"SphericalBox", @"SphereVolume %lf. SphereExtremes %@.",
                         sphereVolume, sphereExtremes);
}

- (void) _useEnvironmentDefaults
{
      seed = 500;
      sphereRadius = 30;
      solventDensity = 1000.0;
      solventDensity *= DENSITY_FACTOR;
}

00465 - (id) initWithEnvironment: (id) object observe: (BOOL) value
{
      if(self = [super initWithEnvironment: object observe: value])
      {
            memoryManager = [AdMemoryManager appMemoryManager];
            dependantsDict = [NSMutableDictionary dictionaryWithCapacity:1];
            [dependantsDict setObject: @"_SphereVolume" forKey: @"SphereRadius"];
            [dependantsDict setObject: @"_SphereVolume" forKey: @"SolventDensity"];
            [dependantsDict retain];

            //\note This will have to be changed in the refactoring

            sphereCentre.vector = 0.0;                            
            sphereCentre.vector = 0.0;
            sphereCentre.vector = 0.0;
            memento = NO;
            currentCaptureMethod = @"Standard";
            
            if(environment != nil)
            {
                  [self synchroniseWithEnvironment];
                  [self registerWithEnvironment];
            }
            else
                  [self _useEnvironmentDefaults];
                                    
            NSDebugLLog(@"SphericalBox", @"SolventDensity %lf. SphereRadius %lf", solventDensity, sphereRadius);
            NSDebugLLog(@"SphericalBox", @"Converted density %lf", solventDensity);
            twister = gsl_rng_alloc(gsl_rng_mt19937);
            gsl_rng_set(twister, seed);
            [self _initialiseDependants];
            [self _createGrid];
      }

      return self;
}

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

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

- (void) dealloc
{
      if(memento)
            [memoryManager freeMatrix: currentSolventCoordinates];
      else
      {
            [memoryManager freeMatrix: solventCoordinates];
            [self _clearCurrentSystem];
            [solventGrid release];
            [dependantsDict release];
            [sphereExtremes release];
            [obscuredIndexes release];
            gsl_rng_free(twister);
      }
      
      [super dealloc];
}

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

Extraction 

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

-(void) hideMoleculesWithIndexes: (NSIndexSet*) indexes
{
      int i, j, k, count;
      int currentAtom;
      
      NSDebugLLog(@"SphericalBox", @"Clearing the current system");

      //free all the information associated with current system
      [self _clearCurrentSystem];
      
      currentNumberOfMolecules = totalNumberOfMolecules - [indexes count];
      currentNumberOfAtoms = currentNumberOfMolecules*atomsPerMolecule;
      currentSolventCoordinates = [memoryManager allocateMatrixWithRows: currentNumberOfAtoms 
                              withColumns: 6]; 
      
      NSDebugLLog(@"SphericalBox", @"There are %d current molecules", currentNumberOfMolecules);

      //copy the unobscured molecule coords to the new matrix

      for(count=0, i=0; i<totalNumberOfMolecules; i++)
            if(![indexes containsIndex: i])     
            {
                  currentAtom = i*atomsPerMolecule;
                  for(k=currentAtom; k<currentAtom + atomsPerMolecule; k++)
                  {
                        for(j=0; j<6; j++)
                              currentSolventCoordinates->matrix[count][j] = solventCoordinates->matrix[k][j];
                        count++;
                  }
            }

      //rebuild the interactions

      NSDebugLLog(@"SphericalBox", @"Rebuilt coordinate matrix");

      [self _setCurrentSystem];
}

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

Public Methods

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

- (BOOL) _checkMolecule: (int) molIndex 
            against: (AdMatrix*) exclusionPoints 
            withExclusionRadius: (double) exclusionRadius
{
      int i,j,k, currentAtom;
      Vector3D seperation;

      for(k=0; k<atomsPerMolecule; k++)
      {
            currentAtom = molIndex*atomsPerMolecule + k;
            for(i=0; i<exclusionPoints->no_rows; i++)
            {
                  for(j=0; j<3; j++)
                        seperation.vector[j] =  solventCoordinates->matrix[currentAtom][j] -
                        exclusionPoints->matrix[i][j];

                  Ad3DVectorLength(&seperation);
                  if(seperation.length < exclusionRadius)
                        return YES;
            }                                         
      }

      return NO;
}

//FIXME: Allocating obscuredIndexes here means this method can only be called once!

- (void) setExclusionPoints: (AdMatrix*) exclusionPoints exclusionRadius: (double) exclusionRadius
{
      int i, j;
      Vector3D seperation;

      NSDebugLLog(@"SphericalBox",
            @"Excluding molecules based on exclusion point matrix with exclusion radius of %lf", 
            exclusionRadius);

      obscuredIndexes = [NSMutableIndexSet new];

      //find what solvent molecules have to been excluded using the completeMatrx
      //then build a new currentSolventCoordinates matrix leaving out the excluded molecules

      for(i=0; i<totalNumberOfMolecules; i++)
            if([self _checkMolecule: i against: exclusionPoints withExclusionRadius: exclusionRadius])
                  [obscuredIndexes addIndex: i];

      NSDebugLLog(@"SphericalBox", @"The following atoms are obscured %@", obscuredIndexes); 
      no_obscured_molecules = [obscuredIndexes count];
      NSDebugLLog(@"SphericalBox", @"There are %d obscured molecules", no_obscured_molecules); 
      [self hideMoleculesWithIndexes: obscuredIndexes];
} 

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

Protocols

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

//DataSourceDelegation

00639 - (void) setDataSource: (id) anObject
{
      if([anObject conformsToProtocol: @protocol(AdSystemDataSource)])
            dataSource = anObject;
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"Object %@ does not conform to AdSystemDataSource protocol", 
                  [anObject description]];
}

00649 - (void) reloadData
{
      NSDebugLLog(@"SphericalBox", @"Loading data from data source");

      if(dataSource != nil)
      {
            [self _cleanUp];
            [self _retrieveCoordinates];
            [self _populateBox];
            [self _setCurrentSystem];
      }
      else
            [NSException raise: NSInternalInconsistencyException
                  format: @"No data source has been set (%@)", [self description]];
}

00665 - (id) dataSource
{
      return dataSource;
}

//Grid Delegate

00672 - (double) cavityVolume
{
      return sphereVolume;
}

00677 - (BOOL) isPointInCavity: (double*) point
{
      int i, j;
      Vector3D seperation;
      
      //return yes if its in the sphere otherwise no

      for(j=0; j<3; j++)
            seperation.vector[j] = point[j] - sphereCentre.vector[j];

      Ad3DVectorLength(&seperation);

      if(seperation.length < sphereRadius)
            return YES;
      
      return NO;
}

00695 - (Vector3D*) cavityCentre
{
      return &sphereCentre;
}

00700 - (NSArray*) cavityExtremes
{
      return sphereExtremes;
}

00705 - (BOOL) allowsTranslation
{
      return YES;
}

00710 - (BOOL) handlesTranslation
{
      return YES;
}

00715 - (BOOL) updateOnTranslation
{
      return NO;
}

//AdSystemDataSource

00722 - (NSValue*) objectValueForCoordinates: (id) sender
{
      return [NSValue valueWithPointer: currentSolventCoordinates];
}

00727 - (NSValue*) objectValueForVelocities: (id) sender
{
      return nil;
}

00732 - (NSValue*) objectValueForAccelerations: (id) sender
{
      return nil;
}

00737 - (NSArray*) objectValueForAtomTypes: (id) sender
{
      return solventAtomTypes;
}

00742 - (NSDictionary*) objectValueForBondedInteractions: (id) object
{
      return solventBondedInteractions;
}

00747 - (NSArray*) objectValueForNonbondedInteractions: (id) object
{
      return solventNonbonded;
}

00752 - (NSDictionary*) objectValueForNonbondedInteractionTypes: (id) object
{
      return solventNonbondedTypes;
}

/*
 * Environment observation
 */

//FIXME: Incomplete implementation of updateForKey:value:object
//Variable dependants not recalulated

00764 - (void) updateForKey: (NSString*) key value: (id) value object: (id) object
{
      if([key isEqual: @"SolvationSphereRadius"])
            sphereRadius = [value doubleValue];
      else if([key isEqual: @"SolventDensity"])
      {
            solventDensity = [value doubleValue];
            solventDensity *= DENSITY_FACTOR;
      }
      else if([key isEqual: @"Seed"])
            seed = [value intValue];
}

00777 - (void) registerWithEnvironment
{
      if(observesEnvironment)
      {     
            [environment addObserver: self forKey: @"SolvationSphereRadius"];
            [environment addObserver: self forKey: @"SolventDensity"];
            [environment addObserver: self forKey: @"Seed"];
      }
}

00787 - (void) deregisterWithEnvironment
{
      [environment removeObserver: self forKey: @"SolvationSphereRadius"];
      [environment removeObserver: self forKey: @"SolventDensity"];
      [environment removeObserver: self forKey: @"Seed"];
}

00794 - (void) synchroniseWithEnvironment
{
      seed = [[environment valueForKey: @"Seed"] intValue];
      sphereRadius =  [[environment valueForKey: @"SolvationSphereRadius"] doubleValue];
      solventDensity = [[environment valueForKey: @"SolventDensity"] doubleValue];
      //FIXME: this should be moved - all units should enter the core converted
      solventDensity *= DENSITY_FACTOR;
}

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

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

Coding

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

- (void) _mementoDecode: (NSCoder*) decoder
{
      int i, j;
      int length, count;
      double* array, *buffer;

      //coordinates

      memento = YES;
      currentSolventCoordinates = (AdMatrix*)malloc(sizeof(AdMatrix));
      buffer = (double*)[decoder decodeBytesForKey: @"CurrentCoordinateMatrix" returnedLength: &length];
      currentSolventCoordinates->no_rows = [decoder decodeIntForKey: @"CurrentCoordinates.Rows"];
      currentSolventCoordinates->no_columns = 6;
      currentSolventCoordinates->matrix = (double**)malloc(currentSolventCoordinates->no_rows*sizeof(double*));
      array = calloc(currentSolventCoordinates->no_rows*6, sizeof(double));
      for(i=0, j=0; i<currentSolventCoordinates->no_rows; i ++, j = j + currentSolventCoordinates->no_columns)
            currentSolventCoordinates->matrix[i] = array + j;

      //The memento only has three columns!

      for(count=0,i=0; i< currentSolventCoordinates->no_rows; i++)
            for(j=0; j<3; j++)
            {
                  currentSolventCoordinates->matrix = buffer;
                  count++;
            }
}

- (void) _fullDecode: (NSCoder*) decoder
{     
      int i, j;
      int length, count;
      double *buffer;
      AdMatrix* dataSourceCoordinates;

      memento = NO;
      dataSource = [decoder decodeObjectForKey: @"DataSource"];
      solventAtomTypes = [[decoder decodeObjectForKey: @"AtomTypes"] retain];

      //coordinates

      buffer = (double*)[decoder decodeBytesForKey: @"CurrentCoordinateMatrix" 
                        returnedLength: &length];
      currentSolventCoordinates = [memoryManager allocateMatrixWithRows: 
                                    [decoder decodeIntForKey: @"CurrentCoordinates.Rows"]
                              withColumns: 6];
      
      for(count=0,i=0; i< currentSolventCoordinates->no_rows; i++)
            for(j=0; j<6; j++)
            {
                  currentSolventCoordinates->matrix = buffer;
                  count++;
            }
      
      buffer = (double*)[decoder decodeBytesForKey: @"CoordinateMatrix" returnedLength: &length];
      solventCoordinates = [memoryManager allocateMatrixWithRows:
                               [decoder decodeIntForKey: @"Coordinates.Rows"]
                        withColumns: 6];

      for(count=0,i=0; i< solventCoordinates->no_rows; i++)
            for(j=0; j<6; j++)
            {
                  solventCoordinates->matrix = buffer;
                  count++;
            }
      
      NSDebugLLog(@"SphericalBox", @"DataSource is %@", dataSource);

      dataSourceCoordinates = [[dataSource objectValueForCoordinates: self] pointerValue];
      for(solventMass=0, i=0; i<dataSourceCoordinates->no_rows; i++)
            solventMass += dataSourceCoordinates->matrix[i][3]; 
      
      //FIXME: If the exact same environment is not present then all these numbers
      //will be meaningless. Objects should only synchronize with those environment variables
      //that it uses dynamically not those it uses for creation etc.

      totalNumberOfMolecules = solventDensity*sphereVolume/solventMass;
      totalNumberOfAtoms = totalNumberOfMolecules*dataSourceCoordinates->no_rows;
      atomsPerMolecule = dataSourceCoordinates->no_rows;
      currentNumberOfAtoms = currentSolventCoordinates->no_rows;
      currentNumberOfMolecules = currentNumberOfAtoms/atomsPerMolecule;
      no_obscured_molecules = totalNumberOfMolecules - currentNumberOfMolecules;

      /*GSPrintf(stderr, @"There are %d molecules, (%d total, %d obscured)\n",
             currentNumberOfMolecules, totalNumberOfMolecules, no_obscured_molecules);*/

      //interactions

      [self _retrieveBondedInteractions];
      [self _retrieveNonbondedInteractions];

      solventGrid = [[decoder decodeObjectForKey: @"SolventGrid"] retain];
      obscuredIndexes = [[decoder decodeObjectForKey: @"ObscuredMolecules"] retain];
      
      //FIXME: Should these be encoded for the memento aswell
      //FIXME: This has to be changed when we introduce moving of the sphere

      sphereCentre.vector = 0.0;                            
      sphereCentre.vector = 0.0;
      sphereCentre.vector = 0.0;
      memento = NO;
      currentCaptureMethod = @"Standard";
      twister = gsl_rng_alloc(gsl_rng_mt19937);
      gsl_rng_set(twister, seed);
}

//FIXME: We dont encode dependantsDict - Have to decide the mechanism
//for update dependant variables

- (id) initWithCoder: (NSCoder*) decoder
{
      NSString* archiveType;

      self = [super initWithCoder: decoder];
      if([decoder allowsKeyedCoding])
      {
            archiveType = [decoder decodeObjectForKey: @"ArchiveType"];
            memoryManager = [AdMemoryManager appMemoryManager];
            environment = [AdEnvironment globalEnvironment];
            if(environment != nil)
            {
                  [self synchroniseWithEnvironment];
                  [self registerWithEnvironment];
            }
            else
                  [self _useEnvironmentDefaults];

            [self _initialiseDependants];

            NSDebugLLog(@"SphericalBox", @"Archive type is %@", archiveType);

            if([archiveType isEqual: @"Memento"])
                  [self _mementoDecode: decoder];
            else
                  [self _fullDecode: decoder];
      }
      else
            [NSException raise: NSInvalidArgumentException 
                  format: @"%@ does not support non keyed coding", [self classDescription]];

      return self;
}

- (void) _mementoEncodeWithCoder: (NSCoder*) encoder
{
      int i, j, count;
      int bytesLength;
      double *buffer;

      NSDebugLLog(@"Encode", @"Encoding %@", [self description]);
      bytesLength = sizeof(AdMatrixSize)*currentSolventCoordinates->no_rows*3;      
      buffer = (double*)malloc(bytesLength*sizeof(double));
      
      for(i=0, count=0; i<currentSolventCoordinates->no_rows; i++)
            for(j=0; j<3; j++)
            {
                  buffer = currentSolventCoordinates->matrix;
                  count++;
            }

      [encoder encodeBytes: (uint8_t*)buffer
            length: bytesLength 
            forKey: @"CurrentCoordinateMatrix"];
      [encoder encodeInt: currentSolventCoordinates->no_rows 
            forKey: @"CurrentCoordinates.Rows"];      
      free(buffer);
}

- (void) _fullEncodeWithCoder: (NSCoder*) encoder
{
      int bytesLength;

      NSDebugLLog(@"Encode", @"Encoding %@", [self description]);
            
      [encoder encodeConditionalObject: dataSource forKey: @"DataSource"];
      [encoder encodeObject: solventAtomTypes forKey: @"AtomTypes"];
            
      //total coordinates

      bytesLength = sizeof(AdMatrixSize)*solventCoordinates->no_rows*solventCoordinates->no_columns;  
      [encoder encodeBytes: (uint8_t*)solventCoordinates->matrix 
            length: bytesLength 
            forKey: @"CoordinateMatrix"];
      [encoder encodeInt: solventCoordinates->no_rows forKey: @"Coordinates.Rows"]; 
      NSDebugLLog(@"Encode", @"Encoding %d bytes. Memento offf", bytesLength);

      //current coordinates

      bytesLength = sizeof(AdMatrixSize)*currentSolventCoordinates->no_rows*currentSolventCoordinates->no_columns;      
      [encoder encodeBytes: (uint8_t*)currentSolventCoordinates->matrix 
            length: bytesLength 
            forKey: @"CurrentCoordinateMatrix"];
      [encoder encodeInt: currentSolventCoordinates->no_rows forKey: @"CurrentCoordinates.Rows"];     
      NSDebugLLog(@"Encode", @"Encoding %d bytes", bytesLength);
      NSDebugLLog(@"Encode", @"Complete %@", [self description]);

      //the grid

      [encoder encodeObject: solventGrid forKey: @"SolventGrid"];

      //obscured atoms

      [encoder encodeObject: obscuredIndexes forKey: @"ObscuredMolecules"];
}

- (void) encodeWithCoder: (NSCoder*) encoder
{
      [super encodeWithCoder: encoder];
      if([encoder allowsKeyedCoding])
      {
            NSDebugLLog(@"Encode", @"Encoding %@", [self description]);
            
            if(memento)
            {
                  [encoder encodeObject: @"Memento" forKey: @"ArchiveType"];
                  NSDebugLLog(@"Encode", @"Memento archive");
                  [self _mementoEncodeWithCoder: encoder];
            }
            else
            {
                  [encoder encodeObject: @"Complete" forKey: @"ArchiveType"];
                  NSDebugLLog(@"Encode", @"Complete archive");
                  [self _fullEncodeWithCoder: encoder];
            }
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"%@ class does not support non keyed coding", [self class]];
}

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

AdMemento

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

01052 - (id) setCaptureMethod: (NSString*) name
{
      if([name isEqual: @"Standard"] || [name isEqual: @"Full"])
      {
            [currentCaptureMethod release];
            currentCaptureMethod = [name retain];
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"Capture Method %@ is not supported by this object"];    
}

01064 - (void) captureStateWithArchiver: (NSCoder*) archiver key: (NSString*) key
{
      memento = YES;
      [archiver encodeObject: self forKey: key];
      memento = NO;
}

01071 - (void) returnToState: (id) stateMemento
{
      [NSException raise: NSInternalInconsistencyException
            format: [NSString stringWithFormat: 
                  @"Warning method %@ not implemented yet", NSStringFromSelector(_cmd)]];
}

01078 - (id) captureState
{
      [NSException raise: NSInternalInconsistencyException
            format: [NSString stringWithFormat: 
                  @"Warning method %@ not implemented yet", NSStringFromSelector(_cmd)]];
}
/****************

Accessors

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

- (NSValue*) coordinates
{
      return [NSValue valueWithPointer: currentSolventCoordinates];
}

- (int) atomsPerMolecule
{
      return atomsPerMolecule;
}

- (int) numberOfAtoms
{
      return totalNumberOfMolecules;
}

- (int) numberOccludedMolecules
{
      return no_obscured_molecules;
}

- (NSString*) systemName
{
      return [dataSource systemName];
}

@end

Generated by  Doxygen 1.6.0   Back to index