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"

namespace oopse {
MoleculeStamp::MoleculeStamp() {
    DefineParameter(Name, "name");
    
    deprecatedKeywords_.insert("nAtoms");
    deprecatedKeywords_.insert("nBonds");
    deprecatedKeywords_.insert("nBends");
    deprecatedKeywords_.insert("nTorsions");
    deprecatedKeywords_.insert("nRigidBodies");
    deprecatedKeywords_.insert("nCutoffGroups");
    
}

MoleculeStamp::~MoleculeStamp() {

}

bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
    bool ret = addIndexSensitiveStamp(atomStamps_, atom);
    if (!ret) {
        std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName()  << " Molecule\n";
    }
    return ret;
    
}

bool MoleculeStamp::addBondStamp( BondStamp* bond) {
    bondStamps_.push_back(bond);
    return true;
}

bool MoleculeStamp::addBendStamp( BendStamp* bend) {
    bendStamps_.push_back(bend);
    return true;
}

bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
    torsionStamps_.push_back(torsion);
    return true;
}

bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
    bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
    if (!ret) {
        std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName()  << " Molecule\n";
    }
    return ret;
}

bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
    cutoffGroupStamps_.push_back(cutoffgroup);
    return true;
}

bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
    return addIndexSensitiveStamp(fragmentStamps_, fragment);
}
    
void MoleculeStamp::validate() {
    DataHolder::validate();

    std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
    if (ai != atomStamps_.end()) {
        std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
    }

     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
     if (ri != rigidBodyStamps_.end()) {
         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
     }
    
    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
    if (fi != fragmentStamps_.end()) {
        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
    }

    //make sure index is not out of range
    int natoms = getNAtoms();
    for(int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        if (bondStamp->getA() >=  natoms && bondStamp->getB() >= natoms) {
            std::cout << "Error in Molecule " << getName() <<  ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
        }
    }
    for(int i = 0; i < getNBends(); ++i) {
        BendStamp* bendStamp = getBendStamp(i);
        std::vector<int> bendAtoms =  bendStamp->getMembers();
        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
        if (j != bendAtoms.end()) {
            std::cout << "Error in Molecule " << getName();
        }

        if (bendAtoms.size() == 2 && !bendStamp->haveGhostVectorSource()) {
            std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing";
        }
    }    
    for(int i = 0; i < getNBends(); ++i) {
        TorsionStamp* torsionStamp = getTorsionStamp(i);
        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
        if (j != torsionAtoms.end()) {
            std::cout << "Error in Molecule " << getName();
        }
    }
    for(int i = 0; i < getNCutoffGroups(); ++i) {
        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
        if (j != cutoffGroupAtoms.end()) {
            std::cout << "Error in Molecule " << getName();
        }
    }
        
    atom2Rigidbody.resize(natoms); 
    // negative number means atom is a free atom, does not belong to rigidbody
    //every element in atom2Rigidbody has unique negative number at the very beginning
    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
        atom2Rigidbody[i] = -1 - i;
    }

    for (int i = 0; i < getNRigidBodies(); ++i) {
        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
        std::vector<int> members = rbStamp->getMembers();
        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
            atom2Rigidbody[*j] = i;                 
        }
    }
    //make sure atoms belong to same rigidbody do not bond to each other
    for(int i = 0; i < getNBonds(); ++i) {
        BondStamp* bondStamp = getBondStamp(i);
        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()])
            std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << "belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
        }
        
    for(int i = 0; i < getNBends(); ++i) {
        BendStamp* bendStamp = getBendStamp(i);
        std::vector<int> bendAtoms =  bendStamp->getMembers();
        std::vector<int> rigidSet(getNRigidBodies(), 0);
        std::vector<int>::iterator j;
        for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
            int rigidbodyIndex = atom2Rigidbody[*j];
            if (rigidbodyIndex >= 0) {
                ++rigidSet[rigidbodyIndex];
                if (rigidSet[rigidbodyIndex] > 1) {
                    std::cout << "Error in Molecule " << getName() << ": ";
                    //std::cout << "atoms of bend " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
                }
            }
        }
    }      
    for(int i = 0; i < getNTorsions(); ++i) {
        TorsionStamp* torsionStamp = getTorsionStamp(i);
        std::vector<int> torsionAtoms =  torsionStamp->getMembers();
        std::vector<int> rigidSet(getNRigidBodies(), 0);
        std::vector<int>::iterator j;
        for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
            int rigidbodyIndex = atom2Rigidbody[*j];
            if (rigidbodyIndex >= 0) {
                ++rigidSet[rigidbodyIndex];
                if (rigidSet[rigidbodyIndex] > 1) {
                    std::cout << "Error in Molecule " << getName() << ": ";
                    //std::cout << "atoms of torsion " <<  << "belong to same rigidbody " << rigidbodyIndex << "\n";                    
                }
            }
        }
    } 


    //fill in bond information into atom
    fillBondInfo();
    findBends();
    findTorsions();

    int nrigidAtoms = 0;
    for (int i = 0; i < getNRigidBodies(); ++i) {
        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
        nrigidAtoms += rbStamp->getNMembers();
    }
    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;

 }

void MoleculeStamp::fillBondInfo() {
}

void MoleculeStamp::findBends() {

}

void MoleculeStamp::findTorsions() {

}

//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->getRigidBodyStamp(i);
    numAtom = rbStamp->getNMembers();
    for(int j = 0; j < numAtom; j++)
      if (rbStamp->getMemberAt(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->getRigidBodyStamp(rb1);
  natomInRb1 =rbStamp1->getNMembers();

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

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

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

  }

  return jointAtomIndexPair;
}

}
Revision:	2469
Committed:	Fri Dec 2 15:38:03 2005 UTC (18 years, 7 months ago) by tim
File size:	12467 byte(s)
Log Message:	End of the Link --> List Return of the Oject-Oriented replace yacc/lex parser with antlr parser
#	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	tim	2469	namespace oopse {
50			MoleculeStamp::MoleculeStamp() {
51			DefineParameter(Name, "name");
52
53			deprecatedKeywords_.insert("nAtoms");
54			deprecatedKeywords_.insert("nBonds");
55			deprecatedKeywords_.insert("nBends");
56			deprecatedKeywords_.insert("nTorsions");
57			deprecatedKeywords_.insert("nRigidBodies");
58			deprecatedKeywords_.insert("nCutoffGroups");
59
60			}
61	gezelter	1490
62	tim	2469	MoleculeStamp::~MoleculeStamp() {
63	gezelter	1490
64	tim	2469	}
65	gezelter	1490
66	tim	2469	bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
67			bool ret = addIndexSensitiveStamp(atomStamps_, atom);
68			if (!ret) {
69			std::cout << "multiple atoms have the same index: " << atom->getIndex() <<" in " << getName() << " Molecule\n";
70			}
71			return ret;
72
73			}
74	gezelter	1490
75	tim	2469	bool MoleculeStamp::addBondStamp( BondStamp* bond) {
76			bondStamps_.push_back(bond);
77			return true;
78	gezelter	1490	}
79
80	tim	2469	bool MoleculeStamp::addBendStamp( BendStamp* bend) {
81			bendStamps_.push_back(bend);
82			return true;
83			}
84	gezelter	1490
85	tim	2469	bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
86			torsionStamps_.push_back(torsion);
87			return true;
88			}
89	gezelter	1490
90	tim	2469	bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
91			bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
92			if (!ret) {
93			std::cout << "multiple rigidbodies have the same index: " << rigidbody->getIndex() <<" in " << getName() << " Molecule\n";
94			}
95			return ret;
96			}
97	gezelter	1490
98	tim	2469	bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
99			cutoffGroupStamps_.push_back(cutoffgroup);
100			return true;
101	gezelter	1490	}
102
103	tim	2469	bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
104			return addIndexSensitiveStamp(fragmentStamps_, fragment);
105	gezelter	1490	}
106
107	tim	2469	void MoleculeStamp::validate() {
108			DataHolder::validate();
109	gezelter	1490
110	tim	2469	std::vector<AtomStamp>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp>(NULL));
111			if (ai != atomStamps_.end()) {
112			std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
113	gezelter	1490	}
114
115	tim	2469	std::vector<RigidBodyStamp>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp>(NULL));
116			if (ri != rigidBodyStamps_.end()) {
117			std::cout << "Error in Molecule " << getName() << ":rigidBody[" << ri - rigidBodyStamps_.begin()<< "] is missing\n";
118			}
119	gezelter	1490
120	tim	2469	std::vector<FragmentStamp>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp>(NULL));
121			if (fi != fragmentStamps_.end()) {
122			std::cout << "Error in Molecule " << getName() << ":fragment[" << fi - fragmentStamps_.begin()<< "] is missing\n";
123	gezelter	1490	}
124
125	tim	2469	//make sure index is not out of range
126			int natoms = getNAtoms();
127			for(int i = 0; i < getNBonds(); ++i) {
128			BondStamp* bondStamp = getBondStamp(i);
129			if (bondStamp->getA() >= natoms && bondStamp->getB() >= natoms) {
130			std::cout << "Error in Molecule " << getName() << ": bond between " << bondStamp->getA() << " and " << bondStamp->getB() << " is invalid\n";
131			}
132	gezelter	1490	}
133	tim	2469	for(int i = 0; i < getNBends(); ++i) {
134			BendStamp* bendStamp = getBendStamp(i);
135			std::vector<int> bendAtoms = bendStamp->getMembers();
136			std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
137			if (j != bendAtoms.end()) {
138			std::cout << "Error in Molecule " << getName();
139			}
140	gezelter	1490
141	tim	2469	if (bendAtoms.size() == 2 && !bendStamp->haveGhostVectorSource()) {
142			std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing";
143			}
144			}
145			for(int i = 0; i < getNBends(); ++i) {
146			TorsionStamp* torsionStamp = getTorsionStamp(i);
147			std::vector<int> torsionAtoms = torsionStamp ->getMembers();
148			std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
149			if (j != torsionAtoms.end()) {
150			std::cout << "Error in Molecule " << getName();
151			}
152	gezelter	1490	}
153	tim	2469	for(int i = 0; i < getNCutoffGroups(); ++i) {
154			CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
155			std::vector<int> cutoffGroupAtoms = cutoffGroupStamp ->getMembers();
156			std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), natoms-1));
157			if (j != cutoffGroupAtoms.end()) {
158			std::cout << "Error in Molecule " << getName();
159			}
160	gezelter	1490	}
161	tim	2469
162			atom2Rigidbody.resize(natoms);
163			// negative number means atom is a free atom, does not belong to rigidbody
164			//every element in atom2Rigidbody has unique negative number at the very beginning
165			for(int i = 0; i < atom2Rigidbody.size(); ++i) {
166			atom2Rigidbody[i] = -1 - i;
167	gezelter	1490	}
168
169	tim	2469	for (int i = 0; i < getNRigidBodies(); ++i) {
170			RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
171			std::vector<int> members = rbStamp->getMembers();
172			for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
173			atom2Rigidbody[*j] = i;
174			}
175	gezelter	1490	}
176	tim	2469	//make sure atoms belong to same rigidbody do not bond to each other
177			for(int i = 0; i < getNBonds(); ++i) {
178			BondStamp* bondStamp = getBondStamp(i);
179			if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()])
180			std::cout << "Error in Molecule " << getName() << ": "<<"bond between " << bondStamp->getA() << " and " << bondStamp->getB() << "belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
181			}
182
183			for(int i = 0; i < getNBends(); ++i) {
184			BendStamp* bendStamp = getBendStamp(i);
185			std::vector<int> bendAtoms = bendStamp->getMembers();
186			std::vector<int> rigidSet(getNRigidBodies(), 0);
187			std::vector<int>::iterator j;
188			for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
189			int rigidbodyIndex = atom2Rigidbody[*j];
190			if (rigidbodyIndex >= 0) {
191			++rigidSet[rigidbodyIndex];
192			if (rigidSet[rigidbodyIndex] > 1) {
193			std::cout << "Error in Molecule " << getName() << ": ";
194			//std::cout << "atoms of bend " << << "belong to same rigidbody " << rigidbodyIndex << "\n";
195			}
196			}
197			}
198			}
199			for(int i = 0; i < getNTorsions(); ++i) {
200			TorsionStamp* torsionStamp = getTorsionStamp(i);
201			std::vector<int> torsionAtoms = torsionStamp->getMembers();
202			std::vector<int> rigidSet(getNRigidBodies(), 0);
203			std::vector<int>::iterator j;
204			for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
205			int rigidbodyIndex = atom2Rigidbody[*j];
206			if (rigidbodyIndex >= 0) {
207			++rigidSet[rigidbodyIndex];
208			if (rigidSet[rigidbodyIndex] > 1) {
209			std::cout << "Error in Molecule " << getName() << ": ";
210			//std::cout << "atoms of torsion " << << "belong to same rigidbody " << rigidbodyIndex << "\n";
211			}
212			}
213			}
214			}
215	gezelter	1490
216
217	tim	2469	//fill in bond information into atom
218			fillBondInfo();
219			findBends();
220			findTorsions();
221	gezelter	1490
222	tim	2469	int nrigidAtoms = 0;
223			for (int i = 0; i < getNRigidBodies(); ++i) {
224			RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
225			nrigidAtoms += rbStamp->getNMembers();
226	gezelter	1490	}
227	tim	2469	nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
228	gezelter	1490
229	tim	2469	}
230	gezelter	1490
231	tim	2469	void MoleculeStamp::fillBondInfo() {
232	gezelter	1490	}
233
234	tim	2469	void MoleculeStamp::findBends() {
235	gezelter	1490
236			}
237
238	tim	2469	void MoleculeStamp::findTorsions() {
239	gezelter	1490
240			}
241
242			//Function Name: isBondInSameRigidBody
243			//Return true is both atoms of the bond belong to the same rigid body, otherwise return false
244			bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
245			int rbA;
246			int rbB;
247			int consAtomA;
248			int consAtomB;
249
250			if (!isAtomInRigidBody(bond->getA(),rbA, consAtomA))
251			return false;
252
253			if(!isAtomInRigidBody(bond->getB(),rbB, consAtomB) )
254			return false;
255
256			if(rbB == rbA)
257			return true;
258			else
259			return false;
260			}
261
262			// Function Name: isAtomInRigidBody
263			//return false if atom does not belong to a rigid body, otherwise return true
264			bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
265			int whichRigidBody;
266			int consAtomIndex;
267
268			return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
269
270			}
271
272			// Function Name: isAtomInRigidBody
273			//return false if atom does not belong to a rigid body otherwise return true and set whichRigidBody
274			//and consAtomIndex
275			//atomIndex : the index of atom in component
276			//whichRigidBody: the index of rigidbody in component
277			//consAtomIndex: the position of joint atom apears in rigidbody's definition
278			bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
279			RigidBodyStamp* rbStamp;
280			int numRb;
281			int numAtom;
282
283			whichRigidBody = -1;
284			consAtomIndex = -1;
285
286			numRb = this->getNRigidBodies();
287
288			for(int i = 0 ; i < numRb; i++){
289	tim	2469	rbStamp = this->getRigidBodyStamp(i);
290	gezelter	1490	numAtom = rbStamp->getNMembers();
291			for(int j = 0; j < numAtom; j++)
292	tim	2469	if (rbStamp->getMemberAt(j) == atomIndex){
293	gezelter	1490	whichRigidBody = i;
294			consAtomIndex = j;
295			return true;
296			}
297			}
298
299			return false;
300
301			}
302
303			//return the position of joint atom apears in rigidbody's definition
304			//for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
305			//actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
306	gezelter	2204	std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
307	gezelter	1490	RigidBodyStamp* rbStamp1;
308			RigidBodyStamp* rbStamp2;
309			int natomInRb1;
310			int natomInRb2;
311			int atomIndex1;
312			int atomIndex2;
313	gezelter	2204	std::vector<std::pair<int, int> > jointAtomIndexPair;
314	gezelter	1490
315	tim	2469	rbStamp1 = this->getRigidBodyStamp(rb1);
316	gezelter	1490	natomInRb1 =rbStamp1->getNMembers();
317
318	tim	2469	rbStamp2 = this->getRigidBodyStamp(rb2);
319	gezelter	1490	natomInRb2 =rbStamp2->getNMembers();
320
321			for(int i = 0; i < natomInRb1; i++){
322	tim	2469	atomIndex1 = rbStamp1->getMemberAt(i);
323	gezelter	1490
324			for(int j= 0; j < natomInRb1; j++){
325	tim	2469	atomIndex2 = rbStamp2->getMemberAt(j);
326	gezelter	1490
327			if(atomIndex1 == atomIndex2){
328	gezelter	1930	jointAtomIndexPair.push_back(std::make_pair(i, j));
329	gezelter	1490	break;
330			}
331
332	tim	2469	}
333	gezelter	1490
334	tim	2469	}
335	gezelter	1490
336			return jointAtomIndexPair;
337			}
338	tim	2469
339			}