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:	1234
Committed:	Fri Jun 4 03:15:31 2004 UTC (21 years, 4 months ago) by tim
File size:	13482 byte(s)
Log Message:	new rattle algorithm is working
#	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	if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
312	rigidBodies[rigidBodyIndex] = the_rigidbody;
313	else {
314	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	" first rigidBody declaration." );
321	}
322
323	return NULL;
324	}
325
326	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	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	short int no_atom, no_rigidbody, no_cutoffgroup;
400
401	if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
402	" was not given.\n" );
403
404	if( !have_atoms ){
405	return strdup( "MoleculeStamp error. Molecule contains no atoms." );
406	}
407
408	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	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	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	n_integrable = n_atoms;
445	for (i = 0; i < n_rigidbodies; i++)
446	n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
447
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	return NULL;
456	}
457
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	}