src/types/MoleculeStamp.cpp

/*
 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
 *
 * The University of Notre Dame grants you ("Licensee") a
 * non-exclusive, royalty free, license to use, modify and
 * redistribute this software in source and binary code form, provided
 * that the following conditions are met:
 *
 * 1. Acknowledgement of the program authors must be made in any
 *    publication of scientific results based in part on use of the
 *    program.  An acceptable form of acknowledgement is citation of
 *    the article in which the program was described (Matthew
 *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
 *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
 *    Parallel Simulation Engine for Molecular Dynamics,"
 *    J. Comput. Chem. 26, pp. 252-271 (2005))
 *
 * 2. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 *
 * 3. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the
 *    distribution.
 *
 * This software is provided "AS IS," without a warranty of any
 * kind. All express or implied conditions, representations and
 * warranties, including any implied warranty of merchantability,
 * fitness for a particular purpose or non-infringement, are hereby
 * excluded.  The University of Notre Dame and its licensors shall not
 * be liable for any damages suffered by licensee as a result of
 * using, modifying or distributing the software or its
 * derivatives. In no event will the University of Notre Dame or its
 * licensors be liable for any lost revenue, profit or data, or for
 * direct, indirect, special, consequential, incidental or punitive
 * damages, however caused and regardless of the theory of liability,
 * arising out of the use of or inability to use software, even if the
 * University of Notre Dame has been advised of the possibility of
 * such damages.
 */
 
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>

#include "types/MoleculeStamp.hpp"

char MoleculeStamp::errMsg[500];

MoleculeStamp::MoleculeStamp(){
  
  n_atoms = 0;
  n_bonds = 0;
  n_bends = 0;
  n_torsions = 0;
  n_rigidbodies = 0;
  n_cutoffgroups = 0;
  n_integrable = 0;

  unhandled = NULL;
  atoms = NULL;
  bonds = NULL;
  bends = NULL;
  torsions = NULL;
  rigidBodies = NULL;
  cutoffGroups = NULL;

  have_name = 0;
  have_atoms = 0;
  have_bonds = 0;
  have_bends = 0;
  have_torsions = 0;
  have_rigidbodies = 0;
  have_cutoffgroups = 0;

}

MoleculeStamp::~MoleculeStamp(){
  int i;
  
  if( unhandled != NULL) delete unhandled;
  
  if( rigidBodies != NULL ) {
    for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
  }
  delete[] rigidBodies;

  if( cutoffGroups != NULL ) {
    for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
  }
  delete[] cutoffGroups;
  
  if( atoms != NULL ){
    for( i=0; i<n_atoms; i++ ) delete atoms[i];
  }
  delete[] atoms;
  
  if( bonds != NULL ){
    for( i=0; i<n_bonds; i++ ) delete bonds[i];
  }
  delete[] bonds;
  
  if( bends != NULL ){
    for( i=0; i<n_bends; i++ ) delete bends[i];
  }
  delete[] bends;
  
  if( torsions != NULL ){
    for( i=0; i<n_torsions; i++ ) delete torsions[i];
  }
  delete[] torsions;
  

}

char* MoleculeStamp::assignString( char* lhs, char* rhs ){
  
  if( !strcmp( lhs, "name" ) ){
    strcpy( name, rhs );
    have_name = 1;
  }
  else{
    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
    else unhandled->add( lhs, rhs );
    have_extras = 1;
  }
  return NULL;
}

char* MoleculeStamp::assignDouble( char* lhs, double rhs ){
  int i;

  if( !strcmp( lhs, "nAtoms" ) ){
    n_atoms = (int)rhs;
    
    if( have_atoms ){
      sprintf( errMsg,
               "MoleculeStamp error, n_atoms already declared"
               " for molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_atoms = 1;
    atoms = new AtomStamp*[n_atoms];
    for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nBonds" ) ){
    n_bonds = (int)rhs;

    if( have_bonds ){
      sprintf( errMsg,  
               "MoleculeStamp error, n_bonds already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_bonds = 1;
    bonds = new BondStamp*[n_bonds];
    for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nBends" ) ){
    n_bends = (int)rhs;

    if( have_bends ){
      sprintf( errMsg,
               "MoleculeStamp error, n_bends already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_bends = 1;
    bends = new BendStamp*[n_bends];
    for( i=0; i<n_bends; i++ ) bends[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nTorsions" ) ){
    n_torsions = (int)rhs;

    if( have_torsions ){
      sprintf( errMsg,
               "MoleculeStamp error, n_torsions already declared for"
               " molecule: %s\n",
               name );
      return strdup( errMsg );
    }
    have_torsions = 1;
    torsions = new TorsionStamp*[n_torsions];
    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
  }

  else if( !strcmp( lhs, "nRigidBodies" ) ){
    n_rigidbodies = (int)rhs;

    if( have_rigidbodies ){
      sprintf( errMsg,
               "MoleculeStamp error, n_rigidbodies already declared for"
               " molecule: %s\n",
               name );
      return strdup( errMsg );
    }
    have_rigidbodies = 1;
    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
  }

  else if( !strcmp( lhs, "nCutoffGroups" ) ){
    n_cutoffgroups = (int)rhs;

    if( have_cutoffgroups ){
      sprintf( errMsg,
               "MoleculeStamp error, n_cutoffgroups already declared for"
               " molecule: %s\n",
               name );
      return strdup( errMsg );
    }
    have_cutoffgroups = 1;
    cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
    for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
  }
  
  else{
    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
    else unhandled->add( lhs, rhs );
    have_extras = 1;
  }
  return NULL;
}

char*  MoleculeStamp::assignInt( char* lhs, int rhs ){
  int i;
  
  if( !strcmp( lhs, "nAtoms" ) ){
    n_atoms = rhs;
    
    if( have_atoms ){
      sprintf( errMsg,
               "MoleculeStamp error, n_atoms already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_atoms = 1;
    atoms = new AtomStamp*[n_atoms];
    for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nBonds" ) ){
    n_bonds = rhs;

    if( have_bonds ){
      sprintf( errMsg,
               "MoleculeStamp error, n_bonds already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_bonds = 1;
    bonds = new BondStamp*[n_bonds];
    for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nBends" ) ){
    n_bends = rhs;

    if( have_bends ){
      sprintf( errMsg,
               "MoleculeStamp error, n_bends already declared for"
               " molecule: %s\n",
               name );
      return strdup( errMsg );
    }
    have_bends = 1;
    bends = new BendStamp*[n_bends];
    for( i=0; i<n_bends; i++ ) bends[i] = NULL;
  }
  
  else if( !strcmp( lhs, "nTorsions" ) ){
    n_torsions = rhs;

    if( have_torsions ){
      sprintf( errMsg,
               "MoleculeStamp error, n_torsions already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_torsions = 1;
    torsions = new TorsionStamp*[n_torsions];
    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
  }

  else if( !strcmp( lhs, "nRigidBodies" ) ){
    n_rigidbodies = rhs;

    if( have_rigidbodies ){
      sprintf( errMsg,
               "MoleculeStamp error, n_rigidbodies already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_rigidbodies = 1;
    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
  }
  else if( !strcmp( lhs, "nCutoffGroups" ) ){
    n_cutoffgroups = rhs;

    if( have_cutoffgroups ){
      sprintf( errMsg,
               "MoleculeStamp error, n_cutoffgroups already declared for"
               " molecule: %s\n",
               name);
      return strdup( errMsg );
    }
    have_cutoffgroups = 1;
    cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
    for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
  }
  else{
    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
    else unhandled->add( lhs, rhs );
    have_extras = 1;
  }
  return NULL;
}

char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
  
  if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
  else {
    if( have_atoms ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range", 
               atomIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nAtoms not given before"
                       " first atom declaration." );
  }

  return NULL;
}

char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody, 
                                   int rigidBodyIndex ){
  
  if( have_rigidbodies && rigidBodyIndex < n_rigidbodies ) 
    rigidBodies[rigidBodyIndex] = the_rigidbody;
  else {
    if( have_rigidbodies ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range", 
               rigidBodyIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nRigidBodies not given before"
                       " first rigidBody declaration." );
  }
  
  return NULL;
}

char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup, 
                                     int cutoffGroupIndex ){
  
  if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups ) 
    cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
  else {
    if( have_cutoffgroups ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range", 
               cutoffGroupIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nCutoffGroups not given before"
                       " first CutoffGroup declaration." );
  }
  
  return NULL;
}

char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
  
  
  if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
  else{
    if( have_bonds ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range", 
               bondIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nBonds not given before"
                       "first bond declaration." );
  }

  return NULL;
}

char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
  
  
  if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
  else{
    if( have_bends ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nBends range", 
               bendIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nBends not given before"
                       "first bend declaration." );
  }

  return NULL;
}

char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
  
  
  if( have_torsions && torsionIndex < n_torsions ) 
    torsions[torsionIndex] = the_torsion;
  else{
    if( have_torsions ){
      sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range", 
               torsionIndex );
      return strdup( errMsg );
    }
    else return strdup("MoleculeStamp error, nTorsions not given before"
                       "first torsion declaration." );
  }

  return NULL;
}

char* MoleculeStamp::checkMe( void ){
  
  int i;
  short int no_atom, no_rigidbody, no_cutoffgroup;

  if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
                                  " was not given.\n" );
  
  if( !have_atoms ){
    return strdup( "MoleculeStamp error. Molecule contains no atoms." );
  }
  
  no_rigidbody = 0;
  for( i=0; i<n_rigidbodies; i++ ){
    if( rigidBodies[i] == NULL ) no_rigidbody = 1;
  }

  if( no_rigidbody ){
    sprintf( errMsg, 
             "MoleculeStamp error. Not all of the RigidBodies were"
             " declared in molecule \"%s\".\n", name );
    return strdup( errMsg );
  }

  no_cutoffgroup = 0;
  for( i=0; i<n_cutoffgroups; i++ ){
    if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
  }

  if( no_cutoffgroup ){
    sprintf( errMsg, 
             "MoleculeStamp error. Not all of the CutoffGroups were"
             " declared in molecule \"%s\".\n", name );
    return strdup( errMsg );
  }
  
  no_atom = 0;
  for( i=0; i<n_atoms; i++ ){
    if( atoms[i] == NULL ) no_atom = 1;
  }

  if( no_atom ){
    sprintf( errMsg, 
             "MoleculeStamp error. Not all of the atoms were"
             " declared in molecule \"%s\".\n", name );
    return strdup( errMsg );
  }

  n_integrable = n_atoms;
  for (i = 0; i < n_rigidbodies; i++) 
    n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
  
  if (n_integrable <= 0 || n_integrable > n_atoms) {
    sprintf( errMsg, 
             "MoleculeStamp error. n_integrable is either <= 0 or"
             " greater than n_atoms in molecule \"%s\".\n", name );
    return strdup( errMsg );
  }

  return NULL;
}  


//Function Name: isBondInSameRigidBody
//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
  int rbA;
  int rbB;
  int consAtomA;
  int consAtomB;

  if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
    return false;

  if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
    return false;

  if(rbB == rbA)
    return true;
  else
    return false;
}

// Function Name: isAtomInRigidBody 
//return false if atom does not belong to a rigid body, otherwise return true 
bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
  int whichRigidBody;
  int consAtomIndex;

  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
   
}

// Function Name: isAtomInRigidBody 
//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody 
//and consAtomIndex
//atomIndex : the index of atom in component
//whichRigidBody: the index of rigidbody in component
//consAtomIndex:  the position of joint atom apears in  rigidbody's definition
bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
  RigidBodyStamp* rbStamp;
  int numRb;
  int numAtom;

  whichRigidBody = -1;
  consAtomIndex = -1;

  numRb = this->getNRigidBodies();
  
  for(int i = 0 ; i < numRb; i++){
    rbStamp = this->getRigidBody(i);
    numAtom = rbStamp->getNMembers();
    for(int j = 0; j < numAtom; j++)
      if (rbStamp->getMember(j) == atomIndex){
        whichRigidBody = i;
        consAtomIndex = j;
        return true;
      }
  }

  return false;
   
}

//return the position of joint atom apears in  rigidbody's definition
//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
  RigidBodyStamp* rbStamp1;
  RigidBodyStamp* rbStamp2;
  int natomInRb1;
  int natomInRb2;
  int atomIndex1;
  int atomIndex2;
  std::vector<std::pair<int, int> > jointAtomIndexPair;
  
  rbStamp1 = this->getRigidBody(rb1);
  natomInRb1 =rbStamp1->getNMembers();

  rbStamp2 = this->getRigidBody(rb2);
  natomInRb2 =rbStamp2->getNMembers();

  for(int i = 0; i < natomInRb1; i++){
    atomIndex1 = rbStamp1->getMember(i);
      
    for(int j= 0; j < natomInRb1; j++){
      atomIndex2 = rbStamp2->getMember(j);

      if(atomIndex1 == atomIndex2){
        jointAtomIndexPair.push_back(std::make_pair(i, j));
        break;
      }
      
    }//end for(j =0)

  }//end for (i = 0)

  return jointAtomIndexPair;
}
Revision:	2204
Committed:	Fri Apr 15 22:04:00 2005 UTC (19 years, 2 months ago) by gezelter
File size:	16010 byte(s)
Log Message:	xemacs has been drafted to perform our indentation services
#	User	Rev	Content
1	gezelter	2204	/*
2	gezelter	1930	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3			*
4			* The University of Notre Dame grants you ("Licensee") a
5			* non-exclusive, royalty free, license to use, modify and
6			* redistribute this software in source and binary code form, provided
7			* that the following conditions are met:
8			*
9			* 1. Acknowledgement of the program authors must be made in any
10			* publication of scientific results based in part on use of the
11			* program. An acceptable form of acknowledgement is citation of
12			* the article in which the program was described (Matthew
13			* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14			* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15			* Parallel Simulation Engine for Molecular Dynamics,"
16			* J. Comput. Chem. 26, pp. 252-271 (2005))
17			*
18			* 2. Redistributions of source code must retain the above copyright
19			* notice, this list of conditions and the following disclaimer.
20			*
21			* 3. Redistributions in binary form must reproduce the above copyright
22			* notice, this list of conditions and the following disclaimer in the
23			* documentation and/or other materials provided with the
24			* distribution.
25			*
26			* This software is provided "AS IS," without a warranty of any
27			* kind. All express or implied conditions, representations and
28			* warranties, including any implied warranty of merchantability,
29			* fitness for a particular purpose or non-infringement, are hereby
30			* excluded. The University of Notre Dame and its licensors shall not
31			* be liable for any damages suffered by licensee as a result of
32			* using, modifying or distributing the software or its
33			* derivatives. In no event will the University of Notre Dame or its
34			* licensors be liable for any lost revenue, profit or data, or for
35			* direct, indirect, special, consequential, incidental or punitive
36			* damages, however caused and regardless of the theory of liability,
37			* arising out of the use of or inability to use software, even if the
38			* University of Notre Dame has been advised of the possibility of
39			* such damages.
40			*/
41
42	gezelter	1490	#include <stdlib.h>
43			#include <stdio.h>
44			#include <string.h>
45			#include <iostream>
46
47	tim	1492	#include "types/MoleculeStamp.hpp"
48	gezelter	1490
49			char MoleculeStamp::errMsg[500];
50
51			MoleculeStamp::MoleculeStamp(){
52
53			n_atoms = 0;
54			n_bonds = 0;
55			n_bends = 0;
56			n_torsions = 0;
57			n_rigidbodies = 0;
58			n_cutoffgroups = 0;
59			n_integrable = 0;
60
61			unhandled = NULL;
62			atoms = NULL;
63			bonds = NULL;
64			bends = NULL;
65			torsions = NULL;
66			rigidBodies = NULL;
67			cutoffGroups = NULL;
68
69			have_name = 0;
70			have_atoms = 0;
71			have_bonds = 0;
72			have_bends = 0;
73			have_torsions = 0;
74			have_rigidbodies = 0;
75			have_cutoffgroups = 0;
76
77			}
78
79			MoleculeStamp::~MoleculeStamp(){
80			int i;
81
82			if( unhandled != NULL) delete unhandled;
83
84			if( rigidBodies != NULL ) {
85			for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
86			}
87			delete[] rigidBodies;
88
89			if( cutoffGroups != NULL ) {
90			for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
91			}
92			delete[] cutoffGroups;
93
94			if( atoms != NULL ){
95			for( i=0; i<n_atoms; i++ ) delete atoms[i];
96			}
97			delete[] atoms;
98
99			if( bonds != NULL ){
100			for( i=0; i<n_bonds; i++ ) delete bonds[i];
101			}
102			delete[] bonds;
103
104			if( bends != NULL ){
105			for( i=0; i<n_bends; i++ ) delete bends[i];
106			}
107			delete[] bends;
108
109			if( torsions != NULL ){
110			for( i=0; i<n_torsions; i++ ) delete torsions[i];
111			}
112			delete[] torsions;
113
114
115
116
117			}
118
119			char* MoleculeStamp::assignString( char* lhs, char* rhs ){
120
121			if( !strcmp( lhs, "name" ) ){
122			strcpy( name, rhs );
123			have_name = 1;
124			}
125			else{
126			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
127			else unhandled->add( lhs, rhs );
128			have_extras = 1;
129			}
130			return NULL;
131			}
132
133			char* MoleculeStamp::assignDouble( char* lhs, double rhs ){
134			int i;
135
136			if( !strcmp( lhs, "nAtoms" ) ){
137			n_atoms = (int)rhs;
138
139			if( have_atoms ){
140			sprintf( errMsg,
141			"MoleculeStamp error, n_atoms already declared"
142			" for molecule: %s\n",
143			name);
144			return strdup( errMsg );
145			}
146			have_atoms = 1;
147			atoms = new AtomStamp*[n_atoms];
148			for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
149			}
150
151			else if( !strcmp( lhs, "nBonds" ) ){
152			n_bonds = (int)rhs;
153
154			if( have_bonds ){
155			sprintf( errMsg,
156			"MoleculeStamp error, n_bonds already declared for"
157			" molecule: %s\n",
158			name);
159			return strdup( errMsg );
160			}
161			have_bonds = 1;
162			bonds = new BondStamp*[n_bonds];
163			for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
164			}
165
166			else if( !strcmp( lhs, "nBends" ) ){
167			n_bends = (int)rhs;
168
169			if( have_bends ){
170			sprintf( errMsg,
171			"MoleculeStamp error, n_bends already declared for"
172			" molecule: %s\n",
173			name);
174			return strdup( errMsg );
175			}
176			have_bends = 1;
177			bends = new BendStamp*[n_bends];
178			for( i=0; i<n_bends; i++ ) bends[i] = NULL;
179			}
180
181			else if( !strcmp( lhs, "nTorsions" ) ){
182			n_torsions = (int)rhs;
183
184			if( have_torsions ){
185			sprintf( errMsg,
186			"MoleculeStamp error, n_torsions already declared for"
187			" molecule: %s\n",
188			name );
189			return strdup( errMsg );
190			}
191			have_torsions = 1;
192			torsions = new TorsionStamp*[n_torsions];
193			for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
194			}
195
196			else if( !strcmp( lhs, "nRigidBodies" ) ){
197			n_rigidbodies = (int)rhs;
198
199			if( have_rigidbodies ){
200			sprintf( errMsg,
201			"MoleculeStamp error, n_rigidbodies already declared for"
202			" molecule: %s\n",
203			name );
204			return strdup( errMsg );
205			}
206			have_rigidbodies = 1;
207			rigidBodies = new RigidBodyStamp*[n_rigidbodies];
208			for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
209			}
210
211			else if( !strcmp( lhs, "nCutoffGroups" ) ){
212			n_cutoffgroups = (int)rhs;
213
214			if( have_cutoffgroups ){
215			sprintf( errMsg,
216			"MoleculeStamp error, n_cutoffgroups already declared for"
217			" molecule: %s\n",
218			name );
219			return strdup( errMsg );
220			}
221			have_cutoffgroups = 1;
222			cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
223			for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
224			}
225
226			else{
227			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
228			else unhandled->add( lhs, rhs );
229			have_extras = 1;
230			}
231			return NULL;
232			}
233
234			char* MoleculeStamp::assignInt( char* lhs, int rhs ){
235			int i;
236
237			if( !strcmp( lhs, "nAtoms" ) ){
238			n_atoms = rhs;
239
240			if( have_atoms ){
241			sprintf( errMsg,
242			"MoleculeStamp error, n_atoms already declared for"
243			" molecule: %s\n",
244			name);
245			return strdup( errMsg );
246			}
247			have_atoms = 1;
248			atoms = new AtomStamp*[n_atoms];
249			for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
250			}
251
252			else if( !strcmp( lhs, "nBonds" ) ){
253			n_bonds = rhs;
254
255			if( have_bonds ){
256			sprintf( errMsg,
257			"MoleculeStamp error, n_bonds already declared for"
258			" molecule: %s\n",
259			name);
260			return strdup( errMsg );
261			}
262			have_bonds = 1;
263			bonds = new BondStamp*[n_bonds];
264			for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
265			}
266
267			else if( !strcmp( lhs, "nBends" ) ){
268			n_bends = rhs;
269
270			if( have_bends ){
271			sprintf( errMsg,
272			"MoleculeStamp error, n_bends already declared for"
273			" molecule: %s\n",
274			name );
275			return strdup( errMsg );
276			}
277			have_bends = 1;
278			bends = new BendStamp*[n_bends];
279			for( i=0; i<n_bends; i++ ) bends[i] = NULL;
280			}
281
282			else if( !strcmp( lhs, "nTorsions" ) ){
283			n_torsions = rhs;
284
285			if( have_torsions ){
286			sprintf( errMsg,
287			"MoleculeStamp error, n_torsions already declared for"
288			" molecule: %s\n",
289			name);
290			return strdup( errMsg );
291			}
292			have_torsions = 1;
293			torsions = new TorsionStamp*[n_torsions];
294			for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
295			}
296
297			else if( !strcmp( lhs, "nRigidBodies" ) ){
298			n_rigidbodies = rhs;
299
300			if( have_rigidbodies ){
301			sprintf( errMsg,
302			"MoleculeStamp error, n_rigidbodies already declared for"
303			" molecule: %s\n",
304			name);
305			return strdup( errMsg );
306			}
307			have_rigidbodies = 1;
308			rigidBodies = new RigidBodyStamp*[n_rigidbodies];
309			for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
310			}
311			else if( !strcmp( lhs, "nCutoffGroups" ) ){
312			n_cutoffgroups = rhs;
313
314			if( have_cutoffgroups ){
315			sprintf( errMsg,
316			"MoleculeStamp error, n_cutoffgroups already declared for"
317			" molecule: %s\n",
318			name);
319			return strdup( errMsg );
320			}
321			have_cutoffgroups = 1;
322			cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
323			for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
324			}
325			else{
326			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
327			else unhandled->add( lhs, rhs );
328			have_extras = 1;
329			}
330			return NULL;
331			}
332
333			char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
334
335			if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
336			else {
337			if( have_atoms ){
338			sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
339			atomIndex );
340			return strdup( errMsg );
341			}
342			else return strdup("MoleculeStamp error, nAtoms not given before"
343			" first atom declaration." );
344			}
345
346			return NULL;
347			}
348
349			char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
350			int rigidBodyIndex ){
351
352			if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
353			rigidBodies[rigidBodyIndex] = the_rigidbody;
354			else {
355			if( have_rigidbodies ){
356			sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
357			rigidBodyIndex );
358			return strdup( errMsg );
359			}
360			else return strdup("MoleculeStamp error, nRigidBodies not given before"
361			" first rigidBody declaration." );
362			}
363
364			return NULL;
365			}
366
367			char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
368			int cutoffGroupIndex ){
369
370			if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
371			cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
372			else {
373			if( have_cutoffgroups ){
374			sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
375			cutoffGroupIndex );
376			return strdup( errMsg );
377			}
378			else return strdup("MoleculeStamp error, nCutoffGroups not given before"
379			" first CutoffGroup declaration." );
380			}
381
382			return NULL;
383			}
384
385			char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
386
387
388			if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
389			else{
390			if( have_bonds ){
391			sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
392			bondIndex );
393			return strdup( errMsg );
394			}
395			else return strdup("MoleculeStamp error, nBonds not given before"
396			"first bond declaration." );
397			}
398
399			return NULL;
400			}
401
402			char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
403
404
405			if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
406			else{
407			if( have_bends ){
408			sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
409			bendIndex );
410			return strdup( errMsg );
411			}
412			else return strdup("MoleculeStamp error, nBends not given before"
413			"first bend declaration." );
414			}
415
416			return NULL;
417			}
418
419			char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
420
421
422			if( have_torsions && torsionIndex < n_torsions )
423			torsions[torsionIndex] = the_torsion;
424			else{
425			if( have_torsions ){
426			sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
427			torsionIndex );
428			return strdup( errMsg );
429			}
430			else return strdup("MoleculeStamp error, nTorsions not given before"
431			"first torsion declaration." );
432			}
433
434			return NULL;
435			}
436
437			char* MoleculeStamp::checkMe( void ){
438
439			int i;
440			short int no_atom, no_rigidbody, no_cutoffgroup;
441
442			if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
443			" was not given.\n" );
444
445			if( !have_atoms ){
446			return strdup( "MoleculeStamp error. Molecule contains no atoms." );
447			}
448
449			no_rigidbody = 0;
450			for( i=0; i<n_rigidbodies; i++ ){
451			if( rigidBodies[i] == NULL ) no_rigidbody = 1;
452			}
453
454			if( no_rigidbody ){
455			sprintf( errMsg,
456			"MoleculeStamp error. Not all of the RigidBodies were"
457			" declared in molecule \"%s\".\n", name );
458			return strdup( errMsg );
459			}
460
461			no_cutoffgroup = 0;
462			for( i=0; i<n_cutoffgroups; i++ ){
463			if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
464			}
465
466			if( no_cutoffgroup ){
467			sprintf( errMsg,
468			"MoleculeStamp error. Not all of the CutoffGroups were"
469			" declared in molecule \"%s\".\n", name );
470			return strdup( errMsg );
471			}
472
473			no_atom = 0;
474			for( i=0; i<n_atoms; i++ ){
475			if( atoms[i] == NULL ) no_atom = 1;
476			}
477
478			if( no_atom ){
479			sprintf( errMsg,
480			"MoleculeStamp error. Not all of the atoms were"
481			" declared in molecule \"%s\".\n", name );
482			return strdup( errMsg );
483			}
484
485			n_integrable = n_atoms;
486			for (i = 0; i < n_rigidbodies; i++)
487			n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
488
489			if (n_integrable <= 0 \|\| n_integrable > n_atoms) {
490			sprintf( errMsg,
491			"MoleculeStamp error. n_integrable is either <= 0 or"
492			" greater than n_atoms in molecule \"%s\".\n", name );
493			return strdup( errMsg );
494			}
495
496			return NULL;
497			}
498
499
500			//Function Name: isBondInSameRigidBody
501			//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
502			bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
503			int rbA;
504			int rbB;
505			int consAtomA;
506			int consAtomB;
507
508			if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
509			return false;
510
511			if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
512			return false;
513
514			if(rbB == rbA)
515			return true;
516			else
517			return false;
518			}
519
520			// Function Name: isAtomInRigidBody
521			//return false if atom does not belong to a rigid body, otherwise return true
522			bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
523			int whichRigidBody;
524			int consAtomIndex;
525
526			return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
527
528			}
529
530			// Function Name: isAtomInRigidBody
531			//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
532			//and consAtomIndex
533			//atomIndex : the index of atom in component
534			//whichRigidBody: the index of rigidbody in component
535			//consAtomIndex: the position of joint atom apears in rigidbody's definition
536			bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
537			RigidBodyStamp* rbStamp;
538			int numRb;
539			int numAtom;
540
541			whichRigidBody = -1;
542			consAtomIndex = -1;
543
544			numRb = this->getNRigidBodies();
545
546			for(int i = 0 ; i < numRb; i++){
547			rbStamp = this->getRigidBody(i);
548			numAtom = rbStamp->getNMembers();
549			for(int j = 0; j < numAtom; j++)
550			if (rbStamp->getMember(j) == atomIndex){
551			whichRigidBody = i;
552			consAtomIndex = j;
553			return true;
554			}
555			}
556
557			return false;
558
559			}
560
561			//return the position of joint atom apears in rigidbody's definition
562			//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
563			//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
564	gezelter	2204	std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
565	gezelter	1490	RigidBodyStamp* rbStamp1;
566			RigidBodyStamp* rbStamp2;
567			int natomInRb1;
568			int natomInRb2;
569			int atomIndex1;
570			int atomIndex2;
571	gezelter	2204	std::vector<std::pair<int, int> > jointAtomIndexPair;
572	gezelter	1490
573			rbStamp1 = this->getRigidBody(rb1);
574			natomInRb1 =rbStamp1->getNMembers();
575
576			rbStamp2 = this->getRigidBody(rb2);
577			natomInRb2 =rbStamp2->getNMembers();
578
579			for(int i = 0; i < natomInRb1; i++){
580			atomIndex1 = rbStamp1->getMember(i);
581
582			for(int j= 0; j < natomInRb1; j++){
583			atomIndex2 = rbStamp2->getMember(j);
584
585			if(atomIndex1 == atomIndex2){
586	gezelter	1930	jointAtomIndexPair.push_back(std::make_pair(i, j));
587	gezelter	1490	break;
588			}
589
590			}//end for(j =0)
591
592			}//end for (i = 0)
593
594			return jointAtomIndexPair;
595			}