OOPSE/libBASS/MoleculeStamp.cpp

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <iostream>

#include "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;

  return isAtomInRigidBody(bond->getA(),rbA, consAtomA) &&
                isAtomInRigidBody(bond->getB(),rbB, consAtomB);
}


// 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;

  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
vector<pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
  RigidBodyStamp* rbStamp1;
  RigidBodyStamp* rbStamp2;
  int natomInRb1;
  int natomInRb2;
  int atomIndex1;
  int atomIndex2;
  vector<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(make_pair(i, j));
        break;
      }
      
    }//end for(j =0)

  }//end for (i = 0)

  return jointAtomIndexPair;
}
Revision:	1257
Committed:	Thu Jun 10 15:00:56 2004 UTC (21 years, 4 months ago) by gezelter
File size:	13482 byte(s)
Log Message:	Removed debugging printf statements
#	User	Rev	Content
1	gezelter	828	#include <stdlib.h>
2			#include <stdio.h>
3			#include <string.h>
4	mmeineke	377	#include <iostream>
5
6			#include "MoleculeStamp.hpp"
7
8			char MoleculeStamp::errMsg[500];
9
10			MoleculeStamp::MoleculeStamp(){
11
12			n_atoms = 0;
13			n_bonds = 0;
14			n_bends = 0;
15			n_torsions = 0;
16	gezelter	957	n_rigidbodies = 0;
17	gezelter	1153	n_cutoffgroups = 0;
18	gezelter	1103	n_integrable = 0;
19	mmeineke	377
20			unhandled = NULL;
21			atoms = NULL;
22			bonds = NULL;
23			bends = NULL;
24			torsions = NULL;
25	gezelter	957	rigidBodies = NULL;
26	gezelter	1153	cutoffGroups = NULL;
27	mmeineke	377
28			have_name = 0;
29			have_atoms = 0;
30			have_bonds = 0;
31			have_bends = 0;
32			have_torsions = 0;
33	gezelter	957	have_rigidbodies = 0;
34	gezelter	1153	have_cutoffgroups = 0;
35	mmeineke	377
36			}
37
38			MoleculeStamp::~MoleculeStamp(){
39			int i;
40
41			if( unhandled != NULL) delete unhandled;
42	gezelter	957
43			if( rigidBodies != NULL ) {
44			for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
45			}
46	gezelter	1216	delete[] rigidBodies;
47	gezelter	1153
48			if( cutoffGroups != NULL ) {
49			for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
50			}
51	gezelter	1216	delete[] cutoffGroups;
52	gezelter	957
53	mmeineke	377	if( atoms != NULL ){
54			for( i=0; i<n_atoms; i++ ) delete atoms[i];
55			}
56	gezelter	1216	delete[] atoms;
57	mmeineke	377
58			if( bonds != NULL ){
59			for( i=0; i<n_bonds; i++ ) delete bonds[i];
60			}
61	gezelter	1216	delete[] bonds;
62	mmeineke	377
63			if( bends != NULL ){
64			for( i=0; i<n_bends; i++ ) delete bends[i];
65			}
66	gezelter	1216	delete[] bends;
67	mmeineke	377
68			if( torsions != NULL ){
69			for( i=0; i<n_torsions; i++ ) delete torsions[i];
70			}
71	gezelter	1216	delete[] torsions;
72	mmeineke	377
73	gezelter	1216
74
75
76	mmeineke	377	}
77
78			char* MoleculeStamp::assignString( char* lhs, char* rhs ){
79
80			if( !strcmp( lhs, "name" ) ){
81			strcpy( name, rhs );
82			have_name = 1;
83			}
84			else{
85			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
86			else unhandled->add( lhs, rhs );
87			have_extras = 1;
88			}
89			return NULL;
90			}
91
92			char* MoleculeStamp::assignDouble( char* lhs, double rhs ){
93			int i;
94
95			if( !strcmp( lhs, "nAtoms" ) ){
96			n_atoms = (int)rhs;
97
98			if( have_atoms ){
99			sprintf( errMsg,
100			"MoleculeStamp error, n_atoms already declared"
101	gezelter	981	" for molecule: %s\n",
102	mmeineke	377	name);
103			return strdup( errMsg );
104			}
105			have_atoms = 1;
106			atoms = new AtomStamp*[n_atoms];
107			for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
108			}
109
110			else if( !strcmp( lhs, "nBonds" ) ){
111			n_bonds = (int)rhs;
112
113			if( have_bonds ){
114			sprintf( errMsg,
115			"MoleculeStamp error, n_bonds already declared for"
116			" molecule: %s\n",
117			name);
118			return strdup( errMsg );
119			}
120			have_bonds = 1;
121			bonds = new BondStamp*[n_bonds];
122			for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
123			}
124
125			else if( !strcmp( lhs, "nBends" ) ){
126			n_bends = (int)rhs;
127
128			if( have_bends ){
129			sprintf( errMsg,
130			"MoleculeStamp error, n_bends already declared for"
131			" molecule: %s\n",
132			name);
133			return strdup( errMsg );
134			}
135			have_bends = 1;
136			bends = new BendStamp*[n_bends];
137			for( i=0; i<n_bends; i++ ) bends[i] = NULL;
138			}
139
140			else if( !strcmp( lhs, "nTorsions" ) ){
141			n_torsions = (int)rhs;
142
143			if( have_torsions ){
144			sprintf( errMsg,
145			"MoleculeStamp error, n_torsions already declared for"
146			" molecule: %s\n",
147			name );
148			return strdup( errMsg );
149			}
150			have_torsions = 1;
151			torsions = new TorsionStamp*[n_torsions];
152			for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
153			}
154	gezelter	957
155			else if( !strcmp( lhs, "nRigidBodies" ) ){
156			n_rigidbodies = (int)rhs;
157
158			if( have_rigidbodies ){
159			sprintf( errMsg,
160			"MoleculeStamp error, n_rigidbodies already declared for"
161			" molecule: %s\n",
162			name );
163			return strdup( errMsg );
164			}
165			have_rigidbodies = 1;
166			rigidBodies = new RigidBodyStamp*[n_rigidbodies];
167			for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
168			}
169
170	gezelter	1153	else if( !strcmp( lhs, "nCutoffGroups" ) ){
171			n_cutoffgroups = (int)rhs;
172
173			if( have_cutoffgroups ){
174			sprintf( errMsg,
175			"MoleculeStamp error, n_cutoffgroups already declared for"
176			" molecule: %s\n",
177			name );
178			return strdup( errMsg );
179			}
180			have_cutoffgroups = 1;
181			cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
182			for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
183			}
184
185	mmeineke	377	else{
186			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
187			else unhandled->add( lhs, rhs );
188			have_extras = 1;
189			}
190			return NULL;
191			}
192
193			char* MoleculeStamp::assignInt( char* lhs, int rhs ){
194			int i;
195
196			if( !strcmp( lhs, "nAtoms" ) ){
197			n_atoms = rhs;
198
199			if( have_atoms ){
200			sprintf( errMsg,
201			"MoleculeStamp error, n_atoms already declared for"
202			" molecule: %s\n",
203			name);
204			return strdup( errMsg );
205			}
206			have_atoms = 1;
207			atoms = new AtomStamp*[n_atoms];
208			for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
209			}
210
211			else if( !strcmp( lhs, "nBonds" ) ){
212			n_bonds = rhs;
213
214			if( have_bonds ){
215			sprintf( errMsg,
216			"MoleculeStamp error, n_bonds already declared for"
217			" molecule: %s\n",
218			name);
219			return strdup( errMsg );
220			}
221			have_bonds = 1;
222			bonds = new BondStamp*[n_bonds];
223			for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
224			}
225
226			else if( !strcmp( lhs, "nBends" ) ){
227			n_bends = rhs;
228
229			if( have_bends ){
230			sprintf( errMsg,
231			"MoleculeStamp error, n_bends already declared for"
232			" molecule: %s\n",
233			name );
234			return strdup( errMsg );
235			}
236			have_bends = 1;
237			bends = new BendStamp*[n_bends];
238			for( i=0; i<n_bends; i++ ) bends[i] = NULL;
239			}
240
241			else if( !strcmp( lhs, "nTorsions" ) ){
242			n_torsions = rhs;
243
244			if( have_torsions ){
245			sprintf( errMsg,
246			"MoleculeStamp error, n_torsions already declared for"
247			" molecule: %s\n",
248			name);
249			return strdup( errMsg );
250			}
251			have_torsions = 1;
252			torsions = new TorsionStamp*[n_torsions];
253			for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
254			}
255	gezelter	957
256			else if( !strcmp( lhs, "nRigidBodies" ) ){
257			n_rigidbodies = rhs;
258
259			if( have_rigidbodies ){
260			sprintf( errMsg,
261	gezelter	1153	"MoleculeStamp error, n_rigidbodies already declared for"
262	gezelter	957	" molecule: %s\n",
263			name);
264			return strdup( errMsg );
265			}
266			have_rigidbodies = 1;
267			rigidBodies = new RigidBodyStamp*[n_rigidbodies];
268			for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
269			}
270	gezelter	1153	else if( !strcmp( lhs, "nCutoffGroups" ) ){
271			n_cutoffgroups = rhs;
272
273			if( have_cutoffgroups ){
274			sprintf( errMsg,
275			"MoleculeStamp error, n_cutoffgroups already declared for"
276			" molecule: %s\n",
277			name);
278			return strdup( errMsg );
279			}
280			have_cutoffgroups = 1;
281			cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
282			for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
283			}
284	mmeineke	377	else{
285			if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
286			else unhandled->add( lhs, rhs );
287			have_extras = 1;
288			}
289			return NULL;
290			}
291
292			char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
293
294			if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
295	gezelter	1103	else {
296	mmeineke	377	if( have_atoms ){
297			sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
298			atomIndex );
299			return strdup( errMsg );
300			}
301			else return strdup("MoleculeStamp error, nAtoms not given before"
302	gezelter	981	" first atom declaration." );
303	mmeineke	377	}
304
305			return NULL;
306			}
307
308	gezelter	957	char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
309			int rigidBodyIndex ){
310
311			if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
312			rigidBodies[rigidBodyIndex] = the_rigidbody;
313	gezelter	1103	else {
314	gezelter	957	if( have_rigidbodies ){
315			sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
316			rigidBodyIndex );
317			return strdup( errMsg );
318			}
319			else return strdup("MoleculeStamp error, nRigidBodies not given before"
320	gezelter	981	" first rigidBody declaration." );
321	gezelter	957	}
322
323			return NULL;
324			}
325
326	gezelter	1153	char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
327			int cutoffGroupIndex ){
328
329			if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
330			cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
331			else {
332			if( have_cutoffgroups ){
333			sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
334			cutoffGroupIndex );
335			return strdup( errMsg );
336			}
337			else return strdup("MoleculeStamp error, nCutoffGroups not given before"
338			" first CutoffGroup declaration." );
339			}
340
341			return NULL;
342			}
343
344	mmeineke	377	char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
345
346
347			if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
348			else{
349			if( have_bonds ){
350			sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
351			bondIndex );
352			return strdup( errMsg );
353			}
354			else return strdup("MoleculeStamp error, nBonds not given before"
355			"first bond declaration." );
356			}
357
358			return NULL;
359			}
360
361			char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
362
363
364			if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
365			else{
366			if( have_bends ){
367			sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
368			bendIndex );
369			return strdup( errMsg );
370			}
371			else return strdup("MoleculeStamp error, nBends not given before"
372			"first bend declaration." );
373			}
374
375			return NULL;
376			}
377
378			char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
379
380
381			if( have_torsions && torsionIndex < n_torsions )
382			torsions[torsionIndex] = the_torsion;
383			else{
384			if( have_torsions ){
385			sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
386			torsionIndex );
387			return strdup( errMsg );
388			}
389			else return strdup("MoleculeStamp error, nTorsions not given before"
390			"first torsion declaration." );
391			}
392
393			return NULL;
394			}
395
396			char* MoleculeStamp::checkMe( void ){
397
398			int i;
399	gezelter	1153	short int no_atom, no_rigidbody, no_cutoffgroup;
400	gezelter	957
401			if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
402			" was not given.\n" );
403	gezelter	988
404			if( !have_atoms ){
405			return strdup( "MoleculeStamp error. Molecule contains no atoms." );
406	mmeineke	377	}
407
408	gezelter	957	no_rigidbody = 0;
409			for( i=0; i<n_rigidbodies; i++ ){
410			if( rigidBodies[i] == NULL ) no_rigidbody = 1;
411			}
412
413			if( no_rigidbody ){
414			sprintf( errMsg,
415			"MoleculeStamp error. Not all of the RigidBodies were"
416			" declared in molecule \"%s\".\n", name );
417			return strdup( errMsg );
418			}
419
420	gezelter	1153	no_cutoffgroup = 0;
421			for( i=0; i<n_cutoffgroups; i++ ){
422			if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
423			}
424
425			if( no_cutoffgroup ){
426			sprintf( errMsg,
427			"MoleculeStamp error. Not all of the CutoffGroups were"
428			" declared in molecule \"%s\".\n", name );
429			return strdup( errMsg );
430			}
431
432	mmeineke	377	no_atom = 0;
433			for( i=0; i<n_atoms; i++ ){
434			if( atoms[i] == NULL ) no_atom = 1;
435			}
436
437			if( no_atom ){
438			sprintf( errMsg,
439			"MoleculeStamp error. Not all of the atoms were"
440			" declared in molecule \"%s\".\n", name );
441			return strdup( errMsg );
442			}
443
444	gezelter	1103	n_integrable = n_atoms;
445			for (i = 0; i < n_rigidbodies; i++)
446	tim	1113	n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
447	gezelter	1103
448			if (n_integrable <= 0 \|\| n_integrable > n_atoms) {
449			sprintf( errMsg,
450			"MoleculeStamp error. n_integrable is either <= 0 or"
451			" greater than n_atoms in molecule \"%s\".\n", name );
452			return strdup( errMsg );
453			}
454
455	mmeineke	377	return NULL;
456	gezelter	988	}
457	tim	1234
458
459			//Function Name: isBondInSameRigidBody
460			//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
461			bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
462			int rbA;
463			int rbB;
464			int consAtomA;
465			int consAtomB;
466
467			return isAtomInRigidBody(bond->getA(),rbA, consAtomA) &&
468			isAtomInRigidBody(bond->getB(),rbB, consAtomB);
469			}
470
471
472			// Function Name isAtomInRigidBody
473			//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
474			//and consAtomIndex
475			//atomIndex : the index of atom in component
476			//whichRigidBody: the index of rigidbody in component
477			//consAtomIndex: the position of joint atom apears in rigidbody's definition
478			bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
479			RigidBodyStamp* rbStamp;
480			int numRb;
481			int numAtom;
482
483			numRb = this->getNRigidBodies();
484
485			for(int i = 0 ; i < numRb; i++){
486			rbStamp = this->getRigidBody(i);
487			numAtom = rbStamp->getNMembers();
488			for(int j = 0; j < numAtom; j++)
489			if (rbStamp->getMember(j) == atomIndex){
490			whichRigidBody = i;
491			consAtomIndex = j;
492			return true;
493			}
494			}
495
496			return false;
497
498			}
499
500			//return the position of joint atom apears in rigidbody's definition
501			//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
502			//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
503			vector<pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
504			RigidBodyStamp* rbStamp1;
505			RigidBodyStamp* rbStamp2;
506			int natomInRb1;
507			int natomInRb2;
508			int atomIndex1;
509			int atomIndex2;
510			vector<pair<int, int> > jointAtomIndexPair;
511
512			rbStamp1 = this->getRigidBody(rb1);
513			natomInRb1 =rbStamp1->getNMembers();
514
515			rbStamp2 = this->getRigidBody(rb2);
516			natomInRb2 =rbStamp2->getNMembers();
517
518			for(int i = 0; i < natomInRb1; i++){
519			atomIndex1 = rbStamp1->getMember(i);
520
521			for(int j= 0; j < natomInRb1; j++){
522			atomIndex2 = rbStamp2->getMember(j);
523
524			if(atomIndex1 == atomIndex2){
525			jointAtomIndexPair.push_back(make_pair(i, j));
526			break;
527			}
528
529			}//end for(j =0)
530
531			}//end for (i = 0)
532
533			return jointAtomIndexPair;
534			}