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

EnzymixForceField.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 "EnzymixForceField.h"

//variables needed by the forcefield functions
//\note to be deprecated

double** accelerations, **coordinates;

@implementation EnzymixForceField

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

Intialisation and clean up methods for various objects
Aswell as retrieval of necessary information from others

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

- (void) _initialiseNonbondedCalculator
{
      NSMutableDictionary *TempDict;

      [nonbondedCalculator setCoordinates: [NSValue valueWithPointer: coordinateMatrix]];
      [nonbondedCalculator setInteractions: [NSValue valueWithPointer: nonbnd_p]];        
      [nonbondedCalculator useExternalForceMatrix: accelerationMatrix];
}

- (void) _initialiseLongRangeNonbondedCalculator
{
      NSWarnLog(@"Method %@ not implemented", NSStringFromSelector(_cmd));
}

- (void) _systemCleanUp
{
      [availableTerms removeAllObjects];
      [[state valueForKey:@"InactiveTerms"] removeAllObjects];

      harmonicBond = harmonicAngle = fourierTorsion = improperTorsion = nonbonded = longRange = NO;
      bnd_pot = ang_pot = tor_pot = vdw_pot = est_pot = itor_pot = total_energy = 0;
}

- (void) _initialisationForSystem
{
      int i, j;

      //FIXME: System Keywords arent defined
      //systemKeywords = [system valueForKey: @"SystemKeywords"];

      //get the coordinates and the accelerations 

      coordinateMatrix = [[system valueForKey: @"coordinates"] pointerValue];
      accelerationMatrix = [[system valueForKey: @"accelerations"] pointerValue];
      no_of_atoms = coordinateMatrix->no_rows;

      //see which of the core terms the system provides information for
      //and check if any longrange nonbonded information is available

      bondedInteractions = [system valueForKey: @"bondedInteractions"];
      nonbondedInteractionTypes = [[system valueForKey: @"nonbondedInteractionTypes"] allKeys];

      if([bondedInteractions objectForKey: @"HarmonicBond"] != nil)
      {
            harmonicBond = YES;
            bonds = (InterTable*)[[bondedInteractions objectForKey: @"HarmonicBond"] pointerValue];
            [availableTerms addObject: @"HarmonicBond"];
      }
      else
            harmonicBond = NO;      
      
      if([bondedInteractions objectForKey: @"HarmonicAngle"] != nil)
      {
            harmonicAngle = YES;
            angles = (InterTable*)[[bondedInteractions objectForKey: @"HarmonicAngle"] pointerValue];
            [availableTerms addObject: @"HarmonicAngle"];
      }
      else
            harmonicAngle = NO;     

      if([bondedInteractions objectForKey: @"FourierTorsion"])
      {
            fourierTorsion = YES;
            torsions = (InterTable*)[[bondedInteractions objectForKey: @"FourierTorsion"] pointerValue];
            [availableTerms addObject: @"FourierTorsion"];
      }     
      else
            fourierTorsion = NO;    
      
      if([bondedInteractions objectForKey: @"HarmonicImproperTorsion"])
      {
            improperTorsion = YES;
            improperTorsions = (InterTable*)[[bondedInteractions 
                  objectForKey: @"HarmonicImproperTorsion"] pointerValue];
            [availableTerms addObject: @"HarmonicImproperTorsion"];
      }     
      else
            improperTorsion = NO;   

      if([nonbondedInteractionTypes containsObject: @"CoulombElectrostatic"] &&
             [nonbondedInteractionTypes containsObject: @"TypeOneVDWInteraction"])
      {
            nonbnd_p = (ListElement*)[[system valueForKey: @"shortRangeNonbondedInteractions"] pointerValue];
            [self _initialiseNonbondedCalculator];
            nonbonded = YES;
            [availableTerms addObject: @"Nonbonded"];
      }     
      else
            nonbonded = NO;   

      //long range nonbonded (we should only retrieve this if the environment indicates .... )
      
      if([nonbondedInteractionTypes containsObject: @"LongRangeNonbonded"])
      {
            longRangeInteractions = (ListElement*)[[system valueForKey: @"longRangeNonbondedInteractions"]
                                           pointerValue];
            [self _initialiseLongRangeNonbondedCalculator];
            [availableTerms addObject: @"LongRangeNonbonded"];
            longRange = YES;
      }     
      else
            longRange = NO;   

      GSPrintf(stderr, @"Available Terms for computation :\n");
      for(i=0; i<(int)[availableTerms count]; i++)
            GSPrintf(stderr, @"\t%@\n", [availableTerms objectAtIndex: i]);
      
      [system setCurrentForceFieldState: state];
}

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

Creation & Maintainence

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

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

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

- (id) initWithEnvironment: (id) object observe: (BOOL) value
{
      NSArray* valueArray, *keyArray;
      NSDictionary *termPotentials;

      if(self = [super initWithEnvironment:object observe: value])
      {
            //we still have to add H-bonds and improper torsions 

            coreTerms = [[NSMutableArray alloc] initWithObjects: 
                              @"HarmonicBond",
                              @"HarmonicAngle",
                              @"FourierTorsion", 
                              @"HarmonicImproperTorsion", 
                              @"TypeOneVDWInteraction",
                              @"CoulombElectrostatic",      
                              nil];

            valueArray = [NSArray arrayWithObjects:
                        [NSValue valueWithPointer: &bnd_pot],
                        [NSValue valueWithPointer: &ang_pot],
                        [NSValue valueWithPointer: &tor_pot],
                        [NSValue valueWithPointer: &itor_pot],
                        [NSValue valueWithPointer: &vdw_pot],
                        [NSValue valueWithPointer: &est_pot],
                        nil];

            termPotentials = [NSDictionary dictionaryWithObjects: valueArray forKeys: coreTerms];           

            keyArray = [NSArray arrayWithObjects: @"ForceField", 
                        @"NonbondedCalculationMethod", 
                        @"LongRangeNonbondedCalculationMethod", 
                        @"InactiveTerms", 
                        @"TermPotentials",
                        @"TotalPotential", 
                        @"CustomTerms", 
                         nil];

            valueArray = [NSArray arrayWithObjects: @"Enzymix", 
                        @"Pure Cutoff",
                        @"",
                        [NSMutableArray arrayWithCapacity: 1],
                        termPotentials, 
                        [NSValue valueWithPointer: &total_energy],
                        [NSMutableDictionary dictionaryWithCapacity: 1],
                        nil];


            customTerms = [[NSMutableDictionary dictionaryWithCapacity: 1] retain];
            state = [[NSMutableDictionary dictionaryWithObjects: valueArray forKeys: keyArray] retain];
            availableTerms = [[NSMutableArray arrayWithCapacity: 1] retain];
            bnd_pot = ang_pot = tor_pot = itor_pot = vdw_pot = est_pot= total_energy = 0;

            if(environment == nil)
            {
                  //create a default nonbondedCalculator 

                  nonbondedCalculator = [PureNonBondedCalculator new];
                  longRangeNonbondedCalculator = nil;
            }
            else
            {
                  nonbondedCalculator = [AdNonbondedCalculator objectForEnvironment: environment];
                  [nonbondedCalculator retain];
                  [state setObject: [environment valueForKey: @"ShortRangeInteractions"]
                         forKey: @"Nonbonded"];
                  
                  //create longRangeNonbondedCalculator

                  longRangeNonbondedCalculator = nil; 

                  //register for observations - nonbonded or longrangeNonbonded calculation method changing
                  
            }
      }

      return self;
}

- (void) dealloc
{
      [state release];
      [availableTerms release];
      [customTerms release];
      [nonbondedCalculator release];
      [super dealloc];
}
      

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

Notifications

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

- (void) handleNotification: (id) SystemNotification
{
      //handle notifications from system regarding bondedInteractions, coordinates or velocities

      NSWarnLog(@"(Enzymix ForceField) Method %@ not implemented", NSStringFromSelector(_cmd));
}


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

Public Methods

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

- (id) calculatePotentialAndUpdateSystem
{
      [self calculatePotential];
      [system setCurrentForceFieldState: state];
      return state;
}

- (id) calculatePotential
{
      register int j;
      double potential;
      NSEnumerator* customTermsEnum;      
      id term, customPotentials;

      bnd_pot = ang_pot = tor_pot = itor_pot = vdw_pot = est_pot = 0;
      coordinates = coordinateMatrix->matrix;
      accelerations = accelerationMatrix->matrix;

      if(nonbonded)
      {
            [nonbondedCalculator evaluatePotential];
            vdw_pot = [nonbondedCalculator VDWPotential];
            est_pot = [nonbondedCalculator ElectrostaticPotential];
      }

      if(harmonicBond)
            for(j=0; j < bonds->no_interactions; j++)
                  AdHarmonicBondEnergy(bonds->table[j], &bnd_pot, coordinates);
      
      if(harmonicAngle)
            for(j=0; j < angles->no_interactions; j++)
                  AdHarmonicAngleEnergy(angles->table[j], &ang_pot, coordinates);
      
      if(fourierTorsion)
            for(j=0; j < torsions->no_interactions; j++)
                  AdFourierTorsionEnergy(torsions->table[j], &tor_pot, coordinates);
      
      if(improperTorsion)
            for(j=0; j < improperTorsions->no_interactions; j++)
                  AdHarmonicImproperTorsionEnergy(improperTorsions->table[j], 
                        &itor_pot, 
                        coordinates);

      if([customTerms count] != 0)
      {
            customPotentials = [state valueForKey: customTerms];
            customTermsEnum = [customTerms keyEnumerator];
            while(term = [customTermsEnum nextObject])
            {     
                  [[customTerms objectForKey: term]  evaluatePotential];
                  potential = [[customTerms objectForKey: term] Potential];   
                  [customPotentials setValue: [NSNumber numberWithDouble: potential]
                        forKey: term];
            }
      }

      total_energy = bnd_pot + ang_pot + tor_pot + vdw_pot + est_pot + itor_pot;
      
      return state;
}


- (void) calculateForces
{
      register int j, i;
      double potential;
      AdMatrix* customForce;
      NSEnumerator* customTermsEnum;      
      id term, customPotentials;
      
      bnd_pot = ang_pot = tor_pot = itor_pot = vdw_pot = est_pot = 0;
      coordinates = coordinateMatrix->matrix;
      accelerations = accelerationMatrix->matrix;

      if(nonbonded)
      {
            [nonbondedCalculator evaluateForces];
            vdw_pot = [nonbondedCalculator VDWPotential];
            est_pot = [nonbondedCalculator ElectrostaticPotential];
      }

      if(harmonicBond)
            for(j=0; j < bonds->no_interactions; j++)
                  AdHarmonicBondForce(bonds->table[j], &bnd_pot);

      if(harmonicAngle)
            for(j=0; j < angles->no_interactions; j++)
                  AdHarmonicAngleForce(angles->table[j], &ang_pot);

      if(fourierTorsion)
            for(j=0; j < torsions->no_interactions; j++)
                  AdFourierTorsionForce(torsions->table[j], &tor_pot);
      
      //Potential Calculation test

      if(improperTorsion)
            for(j=0; j < improperTorsions->no_interactions; j++)
                  AdHarmonicImproperTorsionForce(improperTorsions->table[j], 
                        &itor_pot);

      if([customTerms count] != 0)
      {
            customPotentials = [state valueForKey: customTerms];
            customTermsEnum = [customTerms keyEnumerator];
            while(term = [customTermsEnum nextObject])
            {     
                  [[customTerms objectForKey: term]  evaluateForces];
                  potential = [[customTerms objectForKey: term] Potential];   
                  [customPotentials setValue: [NSNumber numberWithDouble: potential]
                        forKey: term];
                              
                  customForce = [[customTerms objectForKey: term] Forces];
                  for(i=0; i<customForce->no_rows; i++)
                        for(j=0; j<3; j++)
                              accelerations[i][j] += customForce->matrix[i][j];
            }
      }

      total_energy = bnd_pot + ang_pot + tor_pot + vdw_pot + est_pot + itor_pot;
      [system setCurrentForceFieldState: state];
}

//custom term methods

- (void) addCustomTerm: (id) object key: (NSString*) name
{
      if([object isKindOfClass: [AdForceFieldTerm class]])
      {
            [customTerms setObject: object forKey: name];
            [[state valueForKey: @"CustomTerms"] 
                  setObject: [NSNumber numberWithInt: 0]
                  forKey: name];
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"Invalid class for custom force object. Must inherit from AdForceFieldTerm"];
} 

- (void) removeCustomTermForKey: (NSString*) name
{
      [customTerms removeObjectForKey: name];
      [[state valueForKey: @"CustomTerms"]
            removeObjectForKey: name];
}

//Activating/Deactivating Terms

- (void) deactivateTerm: (NSString*) termName
{
      if([availableTerms containsObject: termName])
      {
            if(![[state valueForKey: @"Inactive"] containsObject: termName])
            {
                  if([termName isEqual: @"HarmonicBond"])
                        harmonicBond = NO;
                  if([termName isEqual: @"HarmonicAngle"])
                        harmonicAngle = NO;
                  if([termName isEqual: @"FourierTorsion"])
                        fourierTorsion = NO;
                  if([termName isEqual: @"HarmonicImproperTorsion"])
                        improperTorsion = NO;
                  if([termName isEqual: @"Nonbonded"])
                        nonbonded = NO;
                  if([termName isEqual: @"LongRangeNonbonded"])
                        longRange = NO;

                  [[state valueForKey: @"Inactive"] addObject: termName];
            }
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"Term %@ is not available"];
}

- (void) deactivateTermsWithNames: (NSArray*) names
{
      NSEnumerator* termEnum;
      id term;

      termEnum = [names objectEnumerator];
      while(term = [termEnum nextObject])
            [self deactivateTerm: term];
}

- (void) activateTerm: (NSString*) termName
{
      if([availableTerms containsObject: termName])
      {
            if([[state valueForKey: @"Inactive"] containsObject: termName])
            {
                  if([termName isEqual: @"HarmonicBond"])
                        harmonicBond = YES;
                  if([termName isEqual: @"HarmonicAngle"])
                        harmonicAngle = YES;
                  if([termName isEqual: @"FourierTorsion"])
                        fourierTorsion = YES;
                  if([termName isEqual: @"HarmonicImproperTorsion"])
                        improperTorsion = YES;
                  if([termName isEqual: @"Nonbonded"])
                        nonbonded = YES;
                  if([termName isEqual: @"LongRangeNonbonded"])
                        longRange = YES;

                  [[state valueForKey: @"Inactive"] removeObject: termName];
            }
      }
      else
            [NSException raise: NSInvalidArgumentException
                  format: @"Term %@ is not available"];
}

- (void) activateTermsWithNames: (NSArray*) names
{
      NSEnumerator* termEnum;
      id term;

      termEnum = [names objectEnumerator];
      while(term = [termEnum nextObject])
            [self activateTerm: term];
}


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

Accessors

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

- (id) availableTerms
{
      return [availableTerms copy];
}

- (id) coreTerms
{
      return [coreTerms copy];
}

- (id) system
{
      return system;
}

- (void) setSystem: (id) object
{
      if(system != nil)
      {
            [notificationCenter removeObserver:self name:nil object: system];

            //clean up all variables relating to the last system

            [self _systemCleanUp];
      }

      system = object;

      //register for all notifications from system

      [notificationCenter addObserver: self 
            selector: @selector(handleNotification:)
            name: nil
            object: system];
      
      [self _initialisationForSystem];
}

- (void) setNonbondedCalculator: (id) object
{
      if(![object isKindOfClass: [AdNonbondedCalculator class]])
            [NSException raise: NSInvalidArgumentException
                  format: @"Object does not conform to the AdNonBondedCalculation protocol"];

      [nonbondedCalculator release];
      nonbondedCalculator = object;
      [nonbondedCalculator retain];
      if(system != nil)
            [self _initialiseNonbondedCalculator];
}

- (void) setLongRangeNonbondedCalculator
{
      NSWarnLog(@"Method %@ not implemented", NSStringFromSelector(_cmd));
}

@end

Generated by  Doxygen 1.6.0   Back to index