ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/types/MoleculeStamp.cpp
(Generate patch)

Comparing trunk/OOPSE-2.0/src/types/MoleculeStamp.cpp (file contents):
Revision 1930 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 2507 by tim, Sat Dec 10 16:54:40 2005 UTC

# Line 1 | Line 1
1 < /*
1 > /*
2   * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3   *
4   * The University of Notre Dame grants you ("Licensee") a
# Line 38 | Line 38
38   * University of Notre Dame has been advised of the possibility of
39   * such damages.
40   */
41 –
42 – #include <stdlib.h>
43 – #include <stdio.h>
44 – #include <string.h>
45 – #include <iostream>
41  
42 + #include <functional>
43 + #include <iostream>
44 + #include <sstream>
45   #include "types/MoleculeStamp.hpp"
46 + #include "utils/Tuple.hpp"
47  
48 < char MoleculeStamp::errMsg[500];
48 > namespace oopse {
49  
50 < MoleculeStamp::MoleculeStamp(){
51 <  
52 <  n_atoms = 0;
53 <  n_bonds = 0;
54 <  n_bends = 0;
55 <  n_torsions = 0;
56 <  n_rigidbodies = 0;
57 <  n_cutoffgroups = 0;
58 <  n_integrable = 0;
50 > template<class ContainerType>
51 > std::string containerToString(ContainerType& cont) {
52 >    std::ostringstream oss;
53 >    oss << "(";
54 >    typename ContainerType::iterator i = cont.begin();
55 >    if (i != cont.end()) {
56 >        oss << *i;
57 >        ++i;
58 >    }
59 >    for (; i != cont.end();++i) {
60 >        oss << ", ";
61 >        oss << *i;
62 >    }
63 >    oss << ")";
64 >    return oss.str();
65 > }
66  
67 <  unhandled = NULL;
68 <  atoms = NULL;
69 <  bonds = NULL;
70 <  bends = NULL;
71 <  torsions = NULL;
72 <  rigidBodies = NULL;
73 <  cutoffGroups = NULL;
67 > MoleculeStamp::MoleculeStamp() {
68 >    DefineParameter(Name, "name");
69 >    
70 >    deprecatedKeywords_.insert("nAtoms");
71 >    deprecatedKeywords_.insert("nBonds");
72 >    deprecatedKeywords_.insert("nBends");
73 >    deprecatedKeywords_.insert("nTorsions");
74 >    deprecatedKeywords_.insert("nRigidBodies");
75 >    deprecatedKeywords_.insert("nCutoffGroups");
76 >    
77 > }
78  
79 <  have_name = 0;
70 <  have_atoms = 0;
71 <  have_bonds = 0;
72 <  have_bends = 0;
73 <  have_torsions = 0;
74 <  have_rigidbodies = 0;
75 <  have_cutoffgroups = 0;
79 > MoleculeStamp::~MoleculeStamp() {
80  
81   }
82  
83 < MoleculeStamp::~MoleculeStamp(){
84 <  int i;
85 <  
86 <  if( unhandled != NULL) delete unhandled;
87 <  
88 <  if( rigidBodies != NULL ) {
89 <    for( i=0; i<n_rigidbodies; i++ ) delete rigidBodies[i];
90 <  }
87 <  delete[] rigidBodies;
83 > bool MoleculeStamp::addAtomStamp( AtomStamp* atom) {
84 >    bool ret = addIndexSensitiveStamp(atomStamps_, atom);
85 >    if (!ret) {
86 >         std::cout<< "Error in Molecule " << getName()  << ": multiple atoms have the same indices"<< atom->getIndex() <<"\n";
87 >    }
88 >    return ret;
89 >    
90 > }
91  
92 <  if( cutoffGroups != NULL ) {
93 <    for( i=0; i<n_cutoffgroups; i++ ) delete cutoffGroups[i];
94 <  }
95 <  delete[] cutoffGroups;
93 <  
94 <  if( atoms != NULL ){
95 <    for( i=0; i<n_atoms; i++ ) delete atoms[i];
96 <  }
97 <  delete[] atoms;
98 <  
99 <  if( bonds != NULL ){
100 <    for( i=0; i<n_bonds; i++ ) delete bonds[i];
101 <  }
102 <  delete[] bonds;
103 <  
104 <  if( bends != NULL ){
105 <    for( i=0; i<n_bends; i++ ) delete bends[i];
106 <  }
107 <  delete[] bends;
108 <  
109 <  if( torsions != NULL ){
110 <    for( i=0; i<n_torsions; i++ ) delete torsions[i];
111 <  }
112 <  delete[] torsions;
113 <  
92 > bool MoleculeStamp::addBondStamp( BondStamp* bond) {
93 >    bondStamps_.push_back(bond);
94 >    return true;
95 > }
96  
97 + bool MoleculeStamp::addBendStamp( BendStamp* bend) {
98 +    bendStamps_.push_back(bend);
99 +    return true;
100 + }
101  
102 + bool MoleculeStamp::addTorsionStamp( TorsionStamp* torsion) {
103 +    torsionStamps_.push_back(torsion);
104 +    return true;
105 + }
106  
107 + bool MoleculeStamp::addRigidBodyStamp( RigidBodyStamp* rigidbody) {
108 +    bool ret = addIndexSensitiveStamp(rigidBodyStamps_, rigidbody);
109 +    if (!ret) {
110 +        std::cout<< "Error in Molecule " << getName()  << ": multiple rigidbodies have the same indices: " << rigidbody->getIndex() <<"\n";
111 +    }
112 +    return ret;
113   }
114  
115 < char* MoleculeStamp::assignString( char* lhs, char* rhs ){
116 <  
117 <  if( !strcmp( lhs, "name" ) ){
122 <    strcpy( name, rhs );
123 <    have_name = 1;
124 <  }
125 <  else{
126 <    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
127 <    else unhandled->add( lhs, rhs );
128 <    have_extras = 1;
129 <  }
130 <  return NULL;
115 > bool MoleculeStamp::addCutoffGroupStamp( CutoffGroupStamp* cutoffgroup) {
116 >    cutoffGroupStamps_.push_back(cutoffgroup);
117 >    return true;
118   }
119  
120 < char* MoleculeStamp::assignDouble( char* lhs, double rhs ){
121 <  int i;
122 <
136 <  if( !strcmp( lhs, "nAtoms" ) ){
137 <    n_atoms = (int)rhs;
120 > bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
121 >    return addIndexSensitiveStamp(fragmentStamps_, fragment);
122 > }
123      
124 <    if( have_atoms ){
125 <      sprintf( errMsg,
141 <               "MoleculeStamp error, n_atoms already declared"
142 <               " for molecule: %s\n",
143 <               name);
144 <      return strdup( errMsg );
145 <    }
146 <    have_atoms = 1;
147 <    atoms = new AtomStamp*[n_atoms];
148 <    for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
149 <  }
150 <  
151 <  else if( !strcmp( lhs, "nBonds" ) ){
152 <    n_bonds = (int)rhs;
124 > void MoleculeStamp::validate() {
125 >    DataHolder::validate();
126  
127 <    if( have_bonds ){
128 <      sprintf( errMsg,  
129 <               "MoleculeStamp error, n_bonds already declared for"
130 <               " molecule: %s\n",
131 <               name);
159 <      return strdup( errMsg );
127 >    atom2Rigidbody.resize(getNAtoms());
128 >    // negative number means atom is a free atom, does not belong to rigidbody
129 >    //every element in atom2Rigidbody has unique negative number at the very beginning
130 >    for(int i = 0; i < atom2Rigidbody.size(); ++i) {
131 >        atom2Rigidbody[i] = -1 - i;
132      }
133 <    have_bonds = 1;
134 <    bonds = new BondStamp*[n_bonds];
135 <    for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
136 <  }
137 <  
138 <  else if( !strcmp( lhs, "nBends" ) ){
167 <    n_bends = (int)rhs;
168 <
169 <    if( have_bends ){
170 <      sprintf( errMsg,
171 <               "MoleculeStamp error, n_bends already declared for"
172 <               " molecule: %s\n",
173 <               name);
174 <      return strdup( errMsg );
133 >    for (int i = 0; i < getNRigidBodies(); ++i) {
134 >        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
135 >        std::vector<int> members = rbStamp->getMembers();
136 >        for(std::vector<int>::iterator j = members.begin(); j != members.end(); ++j) {
137 >            atom2Rigidbody[*j] = i;                
138 >        }
139      }
176 –    have_bends = 1;
177 –    bends = new BendStamp*[n_bends];
178 –    for( i=0; i<n_bends; i++ ) bends[i] = NULL;
179 –  }
180 –  
181 –  else if( !strcmp( lhs, "nTorsions" ) ){
182 –    n_torsions = (int)rhs;
140  
141 <    if( have_torsions ){
142 <      sprintf( errMsg,
143 <               "MoleculeStamp error, n_torsions already declared for"
144 <               " molecule: %s\n",
145 <               name );
146 <      return strdup( errMsg );
141 >    checkAtoms();
142 >    checkBonds();
143 >    fillBondInfo();
144 >    checkBends();
145 >    checkTorsions();
146 >    checkRigidBodies();
147 >    checkCutoffGroups();
148 >    checkFragments();
149 >
150 >    int nrigidAtoms = 0;
151 >    for (int i = 0; i < getNRigidBodies(); ++i) {
152 >        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
153 >        nrigidAtoms += rbStamp->getNMembers();
154      }
155 <    have_torsions = 1;
192 <    torsions = new TorsionStamp*[n_torsions];
193 <    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
194 <  }
155 >    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
156  
157 <  else if( !strcmp( lhs, "nRigidBodies" ) ){
197 <    n_rigidbodies = (int)rhs;
157 > }
158  
159 <    if( have_rigidbodies ){
160 <      sprintf( errMsg,
161 <               "MoleculeStamp error, n_rigidbodies already declared for"
162 <               " molecule: %s\n",
203 <               name );
204 <      return strdup( errMsg );
159 > void MoleculeStamp::checkAtoms() {
160 >    std::vector<AtomStamp*>::iterator ai = std::find(atomStamps_.begin(), atomStamps_.end(), static_cast<AtomStamp*>(NULL));
161 >    if (ai != atomStamps_.end()) {
162 >        std::cout << "Error in Molecule " << getName() << ": atom[" << ai - atomStamps_.begin()<< "] is missing\n";
163      }
206 –    have_rigidbodies = 1;
207 –    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
208 –    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
209 –  }
164  
211 –  else if( !strcmp( lhs, "nCutoffGroups" ) ){
212 –    n_cutoffgroups = (int)rhs;
213 –
214 –    if( have_cutoffgroups ){
215 –      sprintf( errMsg,
216 –               "MoleculeStamp error, n_cutoffgroups already declared for"
217 –               " molecule: %s\n",
218 –               name );
219 –      return strdup( errMsg );
220 –    }
221 –    have_cutoffgroups = 1;
222 –    cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
223 –    for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
224 –  }
225 –  
226 –  else{
227 –    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
228 –    else unhandled->add( lhs, rhs );
229 –    have_extras = 1;
230 –  }
231 –  return NULL;
165   }
166  
167 < char*  MoleculeStamp::assignInt( char* lhs, int rhs ){
168 <  int i;
169 <  
170 <  if( !strcmp( lhs, "nAtoms" ) ){
171 <    n_atoms = rhs;
167 > void MoleculeStamp::checkBonds() {
168 >    //make sure index is not out of range
169 >    int natoms = getNAtoms();
170 >    for(int i = 0; i < getNBonds(); ++i) {
171 >        BondStamp* bondStamp = getBondStamp(i);
172 >        if (bondStamp->getA() >=  natoms && bondStamp->getB() >= natoms) {
173 >            std::cout << "Error in Molecule " << getName() <<  ": bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") is invalid\n";
174 >        }
175 >    }
176      
177 <    if( have_atoms ){
178 <      sprintf( errMsg,
179 <               "MoleculeStamp error, n_atoms already declared for"
180 <               " molecule: %s\n",
181 <               name);
182 <      return strdup( errMsg );
177 >    //make sure bonds are unique
178 >    std::set<std::pair<int, int> > allBonds;
179 >    for(int i = 0; i < getNBonds(); ++i) {
180 >        BondStamp* bondStamp= getBondStamp(i);        
181 >        std::pair<int, int> bondPair(bondStamp->getA(), bondStamp->getB());
182 >        //make sure bondTuple.first is always less than or equal to bondTuple.third
183 >        if (bondPair.first > bondPair.second) {
184 >            std::swap(bondPair.first, bondPair.second);
185 >        }
186 >        
187 >        std::set<std::pair<int, int> >::iterator iter = allBonds.find(bondPair);
188 >        if ( iter != allBonds.end()) {
189 >            std::cout << "Error in Molecule " << getName() << ": " << "bond(" <<iter->first << ", "<< iter->second << ")appears multiple times\n";
190 >        } else {
191 >            allBonds.insert(bondPair);
192 >        }
193      }
194 <    have_atoms = 1;
195 <    atoms = new AtomStamp*[n_atoms];
196 <    for( i=0; i<n_atoms; i++ ) atoms[i] = NULL;
197 <  }
198 <  
199 <  else if( !strcmp( lhs, "nBonds" ) ){
200 <    n_bonds = rhs;
254 <
255 <    if( have_bonds ){
256 <      sprintf( errMsg,
257 <               "MoleculeStamp error, n_bonds already declared for"
258 <               " molecule: %s\n",
259 <               name);
260 <      return strdup( errMsg );
194 >    
195 >    //make sure atoms belong to same rigidbody do not bond to each other
196 >    for(int i = 0; i < getNBonds(); ++i) {
197 >        BondStamp* bondStamp = getBondStamp(i);
198 >        if (atom2Rigidbody[bondStamp->getA()] == atom2Rigidbody[bondStamp->getB()]) {
199 >            std::cout << "Error in Molecule " << getName() << ": "<<"bond(" << bondStamp->getA() << ", " << bondStamp->getB() << ") belong to same rigidbody " << atom2Rigidbody[bondStamp->getA()] << "\n";
200 >        }
201      }
202 <    have_bonds = 1;
203 <    bonds = new BondStamp*[n_bonds];
264 <    for( i=0; i<n_bonds; i++ ) bonds[i] = NULL;
265 <  }
266 <  
267 <  else if( !strcmp( lhs, "nBends" ) ){
268 <    n_bends = rhs;
202 >    
203 > }
204  
205 <    if( have_bends ){
206 <      sprintf( errMsg,
207 <               "MoleculeStamp error, n_bends already declared for"
208 <               " molecule: %s\n",
209 <               name );
210 <      return strdup( errMsg );
211 <    }
212 <    have_bends = 1;
278 <    bends = new BendStamp*[n_bends];
279 <    for( i=0; i<n_bends; i++ ) bends[i] = NULL;
280 <  }
281 <  
282 <  else if( !strcmp( lhs, "nTorsions" ) ){
283 <    n_torsions = rhs;
205 > void MoleculeStamp::checkBends() {
206 >    for(int i = 0; i < getNBends(); ++i) {
207 >        BendStamp* bendStamp = getBendStamp(i);
208 >        std::vector<int> bendAtoms =  bendStamp->getMembers();
209 >        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
210 >        if (j != bendAtoms.end()) {
211 >            std::cout << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << "have invalid indices\n";
212 >        }
213  
214 <    if( have_torsions ){
215 <      sprintf( errMsg,
216 <               "MoleculeStamp error, n_torsions already declared for"
217 <               " molecule: %s\n",
218 <               name);
219 <      return strdup( errMsg );
214 >        if (bendAtoms.size() == 2 ) {
215 >            if (!bendStamp->haveGhostVectorSource()) {
216 >                std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
217 >            }else{
218 >                int ghostIndex = bendStamp->getGhostVectorSource();
219 >                if (ghostIndex < getNAtoms()) {
220 >                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
221 >                      std::cout <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
222 >                    }
223 >                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
224 >                        std::cout <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
225 >                    }
226 >                }else {
227 >                    std::cout << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
228 >                }
229 >            }
230 >        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
231 >            std::cout <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
232 >        }
233      }
292 –    have_torsions = 1;
293 –    torsions = new TorsionStamp*[n_torsions];
294 –    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
295 –  }
234  
235 <  else if( !strcmp( lhs, "nRigidBodies" ) ){
236 <    n_rigidbodies = rhs;
237 <
238 <    if( have_rigidbodies ){
239 <      sprintf( errMsg,
240 <               "MoleculeStamp error, n_rigidbodies already declared for"
241 <               " molecule: %s\n",
242 <               name);
243 <      return strdup( errMsg );
244 <    }
245 <    have_rigidbodies = 1;
246 <    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
247 <    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
248 <  }
249 <  else if( !strcmp( lhs, "nCutoffGroups" ) ){
250 <    n_cutoffgroups = rhs;
251 <
252 <    if( have_cutoffgroups ){
253 <      sprintf( errMsg,
254 <               "MoleculeStamp error, n_cutoffgroups already declared for"
255 <               " molecule: %s\n",
256 <               name);
257 <      return strdup( errMsg );
258 <    }
259 <    have_cutoffgroups = 1;
260 <    cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
261 <    for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
262 <  }
263 <  else{
264 <    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
265 <    else unhandled->add( lhs, rhs );
266 <    have_extras = 1;
267 <  }
268 <  return NULL;
269 < }
270 <
271 < char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
272 <  
273 <  if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
274 <  else {
275 <    if( have_atoms ){
276 <      sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
277 <               atomIndex );
278 <      return strdup( errMsg );
279 <    }
280 <    else return strdup("MoleculeStamp error, nAtoms not given before"
281 <                       " first atom declaration." );
282 <  }
235 >    for(int i = 0; i < getNBends(); ++i) {
236 >        BendStamp* bendStamp = getBendStamp(i);
237 >        std::vector<int> bendAtoms =  bendStamp->getMembers();
238 >        std::vector<int> rigidSet(getNRigidBodies(), 0);
239 >        std::vector<int>::iterator j;
240 >        for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
241 >            int rigidbodyIndex = atom2Rigidbody[*j];
242 >            if (rigidbodyIndex >= 0) {
243 >                ++rigidSet[rigidbodyIndex];
244 >                if (rigidSet[rigidbodyIndex] > 1) {
245 >                    std::cout << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";                    
246 >                }
247 >            }
248 >        }
249 >    }
250 >    
251 >    
252 >    std::set<IntTuple3> allBends;
253 >    std::set<IntTuple3>::iterator iter;
254 >    for(int i = 0; i < getNBends(); ++i) {
255 >        BendStamp* bendStamp= getBendStamp(i);
256 >        std::vector<int> bend = bendStamp->getMembers();
257 >        if (bend.size() == 2) {
258 >        // in case we have two ghost bend. For example,
259 >        // bend {
260 >        // members (0, 1);
261 >        //   ghostVectorSource = 0;
262 >        // }
263 >        // and
264 >        // bend {
265 >        //   members (0, 1);
266 >        // ghostVectorSource = 0;
267 >        // }
268 >        // In order to distinguish them. we expand them to Tuple3.
269 >        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
270 >             int ghostIndex = bendStamp->getGhostVectorSource();
271 >             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
272 >             if (j != bend.end()) {
273 >                bend.insert(j, ghostIndex);
274 >             }
275 >        }
276 >        
277 >        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
278 >        //make sure bendTuple.first is always less than or equal to bendTuple.third
279 >        if (bendTuple.first > bendTuple.third) {
280 >            std::swap(bendTuple.first, bendTuple.third);
281 >        }
282 >        
283 >        iter = allBends.find(bendTuple);
284 >        if ( iter != allBends.end()) {
285 >            std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
286 >        } else {
287 >            allBends.insert(bendTuple);
288 >        }
289 >    }
290  
291 <  return NULL;
292 < }
291 >    for (int i = 0; i < getNBonds(); ++i) {
292 >        BondStamp* bondStamp = getBondStamp(i);
293 >        int a = bondStamp->getA();
294 >        int b = bondStamp->getB();
295  
296 < char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
297 <                                   int rigidBodyIndex ){
351 <  
352 <  if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
353 <    rigidBodies[rigidBodyIndex] = the_rigidbody;
354 <  else {
355 <    if( have_rigidbodies ){
356 <      sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
357 <               rigidBodyIndex );
358 <      return strdup( errMsg );
359 <    }
360 <    else return strdup("MoleculeStamp error, nRigidBodies not given before"
361 <                       " first rigidBody declaration." );
362 <  }
363 <  
364 <  return NULL;
365 < }
296 >        AtomStamp* atomA = getAtomStamp(a);
297 >        AtomStamp* atomB = getAtomStamp(b);
298  
299 < char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
300 <                                     int cutoffGroupIndex ){
301 <  
302 <  if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
303 <    cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
304 <  else {
305 <    if( have_cutoffgroups ){
306 <      sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
307 <               cutoffGroupIndex );
308 <      return strdup( errMsg );
309 <    }
378 <    else return strdup("MoleculeStamp error, nCutoffGroups not given before"
379 <                       " first CutoffGroup declaration." );
380 <  }
381 <  
382 <  return NULL;
383 < }
299 >        //find bend c--a--b
300 >        AtomStamp::AtomIter ai;
301 >        for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
302 >        {
303 >            if(b == c)
304 >                continue;          
305 >            
306 >            IntTuple3 newBend(c, a, b);
307 >            if (newBend.first > newBend.third) {
308 >                std::swap(newBend.first, newBend.third);
309 >            }
310  
311 < char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
312 <  
313 <  
314 <  if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
315 <  else{
316 <    if( have_bonds ){
317 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
392 <               bondIndex );
393 <      return strdup( errMsg );
394 <    }
395 <    else return strdup("MoleculeStamp error, nBonds not given before"
396 <                       "first bond declaration." );
397 <  }
311 >            if (allBends.find(newBend) == allBends.end() ) {                
312 >                allBends.insert(newBend);
313 >                BendStamp * newBendStamp = new BendStamp();
314 >                newBendStamp->setMembers(newBend);
315 >                addBendStamp(newBendStamp);
316 >            }
317 >        }        
318  
319 <  return NULL;
320 < }
319 >        //find bend a--b--c
320 >        for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
321 >        {
322 >            if(a == c)
323 >                continue;          
324  
325 < char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
326 <  
327 <  
328 <  if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
329 <  else{
330 <    if( have_bends ){
331 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
332 <               bendIndex );
333 <      return strdup( errMsg );
325 >            IntTuple3 newBend( a, b, c);
326 >            if (newBend.first > newBend.third) {
327 >                std::swap(newBend.first, newBend.third);
328 >            }            
329 >            if (allBends.find(newBend) == allBends.end() ) {                
330 >                allBends.insert(newBend);
331 >                BendStamp * newBendStamp = new BendStamp();
332 >                newBendStamp->setMembers(newBend);
333 >                addBendStamp(newBendStamp);
334 >            }
335 >        }        
336      }
412 –    else return strdup("MoleculeStamp error, nBends not given before"
413 –                       "first bend declaration." );
414 –  }
337  
416 –  return NULL;
338   }
339  
340 < char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
341 <  
342 <  
343 <  if( have_torsions && torsionIndex < n_torsions )
344 <    torsions[torsionIndex] = the_torsion;
345 <  else{
346 <    if( have_torsions ){
347 <      sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
427 <               torsionIndex );
428 <      return strdup( errMsg );
340 > void MoleculeStamp::checkTorsions() {
341 >    for(int i = 0; i < getNTorsions(); ++i) {
342 >        TorsionStamp* torsionStamp = getTorsionStamp(i);
343 >        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
344 >        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
345 >        if (j != torsionAtoms.end()) {
346 >            std::cout << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
347 >        }
348      }
349 <    else return strdup("MoleculeStamp error, nTorsions not given before"
350 <                       "first torsion declaration." );
351 <  }
349 >    
350 >    for(int i = 0; i < getNTorsions(); ++i) {
351 >        TorsionStamp* torsionStamp = getTorsionStamp(i);
352 >        std::vector<int> torsionAtoms =  torsionStamp->getMembers();
353 >        std::vector<int> rigidSet(getNRigidBodies(), 0);
354 >        std::vector<int>::iterator j;
355 >        for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
356 >            int rigidbodyIndex = atom2Rigidbody[*j];
357 >            if (rigidbodyIndex >= 0) {
358 >                ++rigidSet[rigidbodyIndex];
359 >                if (rigidSet[rigidbodyIndex] > 1) {
360 >                    std::cout << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";                  
361 >                }
362 >            }
363 >        }
364 >    }    
365  
366 <  return NULL;
367 < }
366 >    std::set<IntTuple4> allTorsions;
367 >    std::set<IntTuple4>::iterator iter;
368 >     for(int i = 0; i < getNTorsions(); ++i) {
369 >         TorsionStamp* torsionStamp= getTorsionStamp(i);
370 >         std::vector<int> torsion = torsionStamp->getMembers();
371 >         if (torsion.size() == 3) {
372 >             int ghostIndex = torsionStamp->getGhostVectorSource();
373 >             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
374 >             if (j != torsion.end()) {
375 >                torsion.insert(j, ghostIndex);
376 >             }
377 >         }
378  
379 < char* MoleculeStamp::checkMe( void ){
380 <  
381 <  int i;
382 <  short int no_atom, no_rigidbody, no_cutoffgroup;
379 >        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
380 >        if (torsionTuple.first > torsionTuple.fourth) {
381 >            std::swap(torsionTuple.first, torsionTuple.fourth);
382 >            std::swap(torsionTuple.second, torsionTuple.third);                    
383 >        }                
384  
385 <  if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
386 <                                  " was not given.\n" );
387 <  
388 <  if( !have_atoms ){
389 <    return strdup( "MoleculeStamp error. Molecule contains no atoms." );
390 <  }
391 <  
449 <  no_rigidbody = 0;
450 <  for( i=0; i<n_rigidbodies; i++ ){
451 <    if( rigidBodies[i] == NULL ) no_rigidbody = 1;
452 <  }
385 >         iter = allTorsions.find(torsionTuple);
386 >         if ( iter == allTorsions.end()) {
387 >            allTorsions.insert(torsionTuple);
388 >         } else {
389 >            std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
390 >         }
391 >     }
392  
393 <  if( no_rigidbody ){
394 <    sprintf( errMsg,
395 <             "MoleculeStamp error. Not all of the RigidBodies were"
396 <             " declared in molecule \"%s\".\n", name );
458 <    return strdup( errMsg );
459 <  }
393 >    for (int i = 0; i < getNBonds(); ++i) {
394 >        BondStamp* bondStamp = getBondStamp(i);
395 >        int b = bondStamp->getA();
396 >        int c = bondStamp->getB();
397  
398 <  no_cutoffgroup = 0;
399 <  for( i=0; i<n_cutoffgroups; i++ ){
463 <    if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
464 <  }
398 >        AtomStamp* atomB = getAtomStamp(b);
399 >        AtomStamp* atomC = getAtomStamp(c);
400  
401 <  if( no_cutoffgroup ){
402 <    sprintf( errMsg,
468 <             "MoleculeStamp error. Not all of the CutoffGroups were"
469 <             " declared in molecule \"%s\".\n", name );
470 <    return strdup( errMsg );
471 <  }
472 <  
473 <  no_atom = 0;
474 <  for( i=0; i<n_atoms; i++ ){
475 <    if( atoms[i] == NULL ) no_atom = 1;
476 <  }
401 >        AtomStamp::AtomIter ai2;
402 >        AtomStamp::AtomIter ai3;
403  
404 <  if( no_atom ){
405 <    sprintf( errMsg,
406 <             "MoleculeStamp error. Not all of the atoms were"
407 <             " declared in molecule \"%s\".\n", name );
482 <    return strdup( errMsg );
483 <  }
404 >        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
405 >        {
406 >            if(a == c)
407 >                continue;
408  
409 <  n_integrable = n_atoms;
410 <  for (i = 0; i < n_rigidbodies; i++)
411 <    n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
412 <  
413 <  if (n_integrable <= 0 || n_integrable > n_atoms) {
414 <    sprintf( errMsg,
415 <             "MoleculeStamp error. n_integrable is either <= 0 or"
416 <             " greater than n_atoms in molecule \"%s\".\n", name );
417 <    return strdup( errMsg );
418 <  }
409 >            for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
410 >            {
411 >                if(d == b)
412 >                    continue;
413 >                
414 >                IntTuple4 newTorsion(a, b, c, d);
415 >                //make sure the first element is always less than or equal to the fourth element in IntTuple4
416 >                if (newTorsion.first > newTorsion.fourth) {
417 >                    std::swap(newTorsion.first, newTorsion.fourth);
418 >                    std::swap(newTorsion.second, newTorsion.third);                    
419 >                }                
420 >                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
421 >                    allTorsions.insert(newTorsion);
422 >                    TorsionStamp * newTorsionStamp = new TorsionStamp();
423 >                    newTorsionStamp->setMembers(newTorsion);
424 >                    addTorsionStamp(newTorsionStamp);                    
425 >                }            
426 >            }
427 >        }    
428 >    }
429  
430 <  return NULL;
497 < }  
430 > }
431  
432 + void MoleculeStamp::checkRigidBodies() {
433 +     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
434 +     if (ri != rigidBodyStamps_.end()) {
435 +         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
436 +     }
437  
438 +    for (int i = 0; i < getNRigidBodies(); ++i) {
439 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
440 +        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
441 +        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
442 +        if (j != rigidAtoms.end()) {
443 +            std::cout << "Error in Molecule " << getName();
444 +        }
445 +        
446 +    }    
447 + }
448 +
449 + void MoleculeStamp::checkCutoffGroups() {
450 +
451 +    for(int i = 0; i < getNCutoffGroups(); ++i) {
452 +        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
453 +        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
454 +        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
455 +        if (j != cutoffGroupAtoms.end()) {
456 +            std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
457 +        }
458 +    }    
459 + }
460 +
461 + void MoleculeStamp::checkFragments() {
462 +
463 +    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
464 +    if (fi != fragmentStamps_.end()) {
465 +        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
466 +    }
467 +    
468 + }
469 +
470 + void MoleculeStamp::fillBondInfo() {
471 +
472 +    for (int i = 0; i < getNBonds(); ++i) {
473 +        BondStamp* bondStamp = getBondStamp(i);
474 +        int a = bondStamp->getA();
475 +        int b = bondStamp->getB();
476 +        AtomStamp* atomA = getAtomStamp(a);
477 +        AtomStamp* atomB = getAtomStamp(b);
478 +        atomA->addBond(i);
479 +        atomA->addBondedAtom(b);
480 +        atomB->addBond(i);        
481 +        atomB->addBondedAtom(a);
482 +
483 +    }
484 + }
485 +
486 +
487 +
488   //Function Name: isBondInSameRigidBody
489   //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
490   bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
# Line 520 | Line 508 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
508   // Function Name: isAtomInRigidBody
509   //return false if atom does not belong to a rigid body, otherwise return true
510   bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
511 <  int whichRigidBody;
524 <  int consAtomIndex;
525 <
526 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
511 >  return atom2Rigidbody[atomIndex] >=0 ;
512    
513   }
514  
# Line 534 | Line 519 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
519   //whichRigidBody: the index of rigidbody in component
520   //consAtomIndex:  the position of joint atom apears in  rigidbody's definition
521   bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
537 –  RigidBodyStamp* rbStamp;
538 –  int numRb;
539 –  int numAtom;
522  
523 +  
524 +
525    whichRigidBody = -1;
526    consAtomIndex = -1;
527  
528 <  numRb = this->getNRigidBodies();
529 <  
530 <  for(int i = 0 ; i < numRb; i++){
531 <    rbStamp = this->getRigidBody(i);
532 <    numAtom = rbStamp->getNMembers();
533 <    for(int j = 0; j < numAtom; j++)
550 <      if (rbStamp->getMember(j) == atomIndex){
551 <        whichRigidBody = i;
528 >  if (atom2Rigidbody[atomIndex] >=0) {
529 >    whichRigidBody = atom2Rigidbody[atomIndex];
530 >    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
531 >    int numAtom = rbStamp->getNMembers();
532 >    for(int j = 0; j < numAtom; j++) {
533 >      if (rbStamp->getMemberAt(j) == atomIndex){
534          consAtomIndex = j;
535          return true;
536        }
537 +    }
538    }
539  
540    return false;
# Line 561 | Line 544 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
544   //return the position of joint atom apears in  rigidbody's definition
545   //for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
546   //actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
547 < std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
547 > std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
548    RigidBodyStamp* rbStamp1;
549    RigidBodyStamp* rbStamp2;
550    int natomInRb1;
551    int natomInRb2;
552    int atomIndex1;
553    int atomIndex2;
554 <   std::vector<std::pair<int, int> > jointAtomIndexPair;
554 >  std::vector<std::pair<int, int> > jointAtomIndexPair;
555    
556 <  rbStamp1 = this->getRigidBody(rb1);
556 >  rbStamp1 = this->getRigidBodyStamp(rb1);
557    natomInRb1 =rbStamp1->getNMembers();
558  
559 <  rbStamp2 = this->getRigidBody(rb2);
559 >  rbStamp2 = this->getRigidBodyStamp(rb2);
560    natomInRb2 =rbStamp2->getNMembers();
561  
562    for(int i = 0; i < natomInRb1; i++){
563 <    atomIndex1 = rbStamp1->getMember(i);
563 >    atomIndex1 = rbStamp1->getMemberAt(i);
564        
565      for(int j= 0; j < natomInRb1; j++){
566 <      atomIndex2 = rbStamp2->getMember(j);
566 >      atomIndex2 = rbStamp2->getMemberAt(j);
567  
568        if(atomIndex1 == atomIndex2){
569          jointAtomIndexPair.push_back(std::make_pair(i, j));
570          break;
571        }
572        
573 <    }//end for(j =0)
573 >    }
574  
575 <  }//end for (i = 0)
575 >  }
576  
577    return jointAtomIndexPair;
578   }
579 +
580 + }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines