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 ){
  

  printf("rigidBodyIndex = %d\n", 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:	1256
Committed:	Thu Jun 10 14:59:14 2004 UTC (20 years, 1 month ago) by gezelter
File size:	13534 byte(s)
Log Message:	Fixed indexing bug in stamps
#	Content
1	#include <stdlib.h>
2	#include <stdio.h>
3	#include <string.h>
4	#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	n_rigidbodies = 0;
17	n_cutoffgroups = 0;
18	n_integrable = 0;
19
20	unhandled = NULL;
21	atoms = NULL;
22	bonds = NULL;
23	bends = NULL;
24	torsions = NULL;
25	rigidBodies = NULL;
26	cutoffGroups = NULL;
27
28	have_name = 0;
29	have_atoms = 0;
30	have_bonds = 0;
31	have_bends = 0;
32	have_torsions = 0;
33	have_rigidbodies = 0;
34	have_cutoffgroups = 0;
35
36	}
37
38	MoleculeStamp::~MoleculeStamp(){
39	int i;
40
41	if( unhandled != NULL) delete unhandled;
42
43	if( rigidBodies != NULL ) {
44	for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
45	}
46	delete[] rigidBodies;
47
48	if( cutoffGroups != NULL ) {
49	for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
50	}
51	delete[] cutoffGroups;
52
53	if( atoms != NULL ){
54	for( i=0; i<n_atoms; i++ ) delete atoms[i];
55	}
56	delete[] atoms;
57
58	if( bonds != NULL ){
59	for( i=0; i<n_bonds; i++ ) delete bonds[i];
60	}
61	delete[] bonds;
62
63	if( bends != NULL ){
64	for( i=0; i<n_bends; i++ ) delete bends[i];
65	}
66	delete[] bends;
67
68	if( torsions != NULL ){
69	for( i=0; i<n_torsions; i++ ) delete torsions[i];
70	}
71	delete[] torsions;
72
73
74
75
76	}
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	" for molecule: %s\n",
102	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
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	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	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
256	else if( !strcmp( lhs, "nRigidBodies" ) ){
257	n_rigidbodies = rhs;
258
259	if( have_rigidbodies ){
260	sprintf( errMsg,
261	"MoleculeStamp error, n_rigidbodies already declared for"
262	" 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	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	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	else {
296	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	" first atom declaration." );
303	}
304
305	return NULL;
306	}
307
308	char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
309	int rigidBodyIndex ){
310
311
312	printf("rigidBodyIndex = %d\n", rigidBodyIndex);
313	if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
314	rigidBodies[rigidBodyIndex] = the_rigidbody;
315	else {
316	if( have_rigidbodies ){
317	sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
318	rigidBodyIndex );
319	return strdup( errMsg );
320	}
321	else return strdup("MoleculeStamp error, nRigidBodies not given before"
322	" first rigidBody declaration." );
323	}
324
325	return NULL;
326	}
327
328	char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
329	int cutoffGroupIndex ){
330
331	if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
332	cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
333	else {
334	if( have_cutoffgroups ){
335	sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
336	cutoffGroupIndex );
337	return strdup( errMsg );
338	}
339	else return strdup("MoleculeStamp error, nCutoffGroups not given before"
340	" first CutoffGroup declaration." );
341	}
342
343	return NULL;
344	}
345
346	char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
347
348
349	if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
350	else{
351	if( have_bonds ){
352	sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
353	bondIndex );
354	return strdup( errMsg );
355	}
356	else return strdup("MoleculeStamp error, nBonds not given before"
357	"first bond declaration." );
358	}
359
360	return NULL;
361	}
362
363	char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
364
365
366	if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
367	else{
368	if( have_bends ){
369	sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
370	bendIndex );
371	return strdup( errMsg );
372	}
373	else return strdup("MoleculeStamp error, nBends not given before"
374	"first bend declaration." );
375	}
376
377	return NULL;
378	}
379
380	char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
381
382
383	if( have_torsions && torsionIndex < n_torsions )
384	torsions[torsionIndex] = the_torsion;
385	else{
386	if( have_torsions ){
387	sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
388	torsionIndex );
389	return strdup( errMsg );
390	}
391	else return strdup("MoleculeStamp error, nTorsions not given before"
392	"first torsion declaration." );
393	}
394
395	return NULL;
396	}
397
398	char* MoleculeStamp::checkMe( void ){
399
400	int i;
401	short int no_atom, no_rigidbody, no_cutoffgroup;
402
403	if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
404	" was not given.\n" );
405
406	if( !have_atoms ){
407	return strdup( "MoleculeStamp error. Molecule contains no atoms." );
408	}
409
410	no_rigidbody = 0;
411	for( i=0; i<n_rigidbodies; i++ ){
412	if( rigidBodies[i] == NULL ) no_rigidbody = 1;
413	}
414
415	if( no_rigidbody ){
416	sprintf( errMsg,
417	"MoleculeStamp error. Not all of the RigidBodies were"
418	" declared in molecule \"%s\".\n", name );
419	return strdup( errMsg );
420	}
421
422	no_cutoffgroup = 0;
423	for( i=0; i<n_cutoffgroups; i++ ){
424	if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
425	}
426
427	if( no_cutoffgroup ){
428	sprintf( errMsg,
429	"MoleculeStamp error. Not all of the CutoffGroups were"
430	" declared in molecule \"%s\".\n", name );
431	return strdup( errMsg );
432	}
433
434	no_atom = 0;
435	for( i=0; i<n_atoms; i++ ){
436	if( atoms[i] == NULL ) no_atom = 1;
437	}
438
439	if( no_atom ){
440	sprintf( errMsg,
441	"MoleculeStamp error. Not all of the atoms were"
442	" declared in molecule \"%s\".\n", name );
443	return strdup( errMsg );
444	}
445
446	n_integrable = n_atoms;
447	for (i = 0; i < n_rigidbodies; i++)
448	n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
449
450	if (n_integrable <= 0 \|\| n_integrable > n_atoms) {
451	sprintf( errMsg,
452	"MoleculeStamp error. n_integrable is either <= 0 or"
453	" greater than n_atoms in molecule \"%s\".\n", name );
454	return strdup( errMsg );
455	}
456
457	return NULL;
458	}
459
460
461	//Function Name: isBondInSameRigidBody
462	//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
463	bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
464	int rbA;
465	int rbB;
466	int consAtomA;
467	int consAtomB;
468
469	return isAtomInRigidBody(bond->getA(),rbA, consAtomA) &&
470	isAtomInRigidBody(bond->getB(),rbB, consAtomB);
471	}
472
473
474	// Function Name isAtomInRigidBody
475	//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
476	//and consAtomIndex
477	//atomIndex : the index of atom in component
478	//whichRigidBody: the index of rigidbody in component
479	//consAtomIndex: the position of joint atom apears in rigidbody's definition
480	bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
481	RigidBodyStamp* rbStamp;
482	int numRb;
483	int numAtom;
484
485	numRb = this->getNRigidBodies();
486
487	for(int i = 0 ; i < numRb; i++){
488	rbStamp = this->getRigidBody(i);
489	numAtom = rbStamp->getNMembers();
490	for(int j = 0; j < numAtom; j++)
491	if (rbStamp->getMember(j) == atomIndex){
492	whichRigidBody = i;
493	consAtomIndex = j;
494	return true;
495	}
496	}
497
498	return false;
499
500	}
501
502	//return the position of joint atom apears in rigidbody's definition
503	//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
504	//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
505	vector<pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
506	RigidBodyStamp* rbStamp1;
507	RigidBodyStamp* rbStamp2;
508	int natomInRb1;
509	int natomInRb2;
510	int atomIndex1;
511	int atomIndex2;
512	vector<pair<int, int> > jointAtomIndexPair;
513
514	rbStamp1 = this->getRigidBody(rb1);
515	natomInRb1 =rbStamp1->getNMembers();
516
517	rbStamp2 = this->getRigidBody(rb2);
518	natomInRb2 =rbStamp2->getNMembers();
519
520	for(int i = 0; i < natomInRb1; i++){
521	atomIndex1 = rbStamp1->getMember(i);
522
523	for(int j= 0; j < natomInRb1; j++){
524	atomIndex2 = rbStamp2->getMember(j);
525
526	if(atomIndex1 == atomIndex2){
527	jointAtomIndexPair.push_back(make_pair(i, j));
528	break;
529	}
530
531	}//end for(j =0)
532
533	}//end for (i = 0)
534
535	return jointAtomIndexPair;
536	}