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

Comparing trunk/OOPSE-3.0/src/types/MoleculeStamp.cpp (file contents):
Revision 1492 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 2493 by tim, Tue Dec 6 21:01:08 2005 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
2 < #include <stdio.h>
3 < #include <string.h>
4 < #include <iostream>
1 > /*
2 > * 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 + #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;
29 <  have_atoms = 0;
30 <  have_bonds = 0;
31 <  have_bends = 0;
32 <  have_torsions = 0;
33 <  have_rigidbodies = 0;
34 <  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 <  }
46 <  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;
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 <  
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" ) ){
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;
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 <
95 <  if( !strcmp( lhs, "nAtoms" ) ){
96 <    n_atoms = (int)rhs;
120 > bool MoleculeStamp::addFragmentStamp( FragmentStamp* fragment) {
121 >    return addIndexSensitiveStamp(fragmentStamps_, fragment);
122 > }
123      
124 <    if( have_atoms ){
125 <      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;
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);
118 <      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" ) ){
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 );
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      }
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;
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 );
147 <    }
148 <    have_torsions = 1;
151 <    torsions = new TorsionStamp*[n_torsions];
152 <    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
153 <  }
141 >    checkAtoms();
142 >    checkBonds();
143 >    fillBondInfo();
144 >    checkBends();
145 >    checkTorsions();
146 >    checkRigidBodies();
147 >    checkCutoffGroups();
148 >    checkFragments();
149  
150 <  else if( !strcmp( lhs, "nRigidBodies" ) ){
151 <    n_rigidbodies = (int)rhs;
152 <
153 <    if( have_rigidbodies ){
159 <      sprintf( errMsg,
160 <               "MoleculeStamp error, n_rigidbodies already declared for"
161 <               " molecule: %s\n",
162 <               name );
163 <      return strdup( errMsg );
150 >    int nrigidAtoms = 0;
151 >    for (int i = 0; i < getNRigidBodies(); ++i) {
152 >        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
153 >        nrigidAtoms += rbStamp->getNMembers();
154      }
155 <    have_rigidbodies = 1;
166 <    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
167 <    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
168 <  }
155 >    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
156  
157 <  else if( !strcmp( lhs, "nCutoffGroups" ) ){
171 <    n_cutoffgroups = (int)rhs;
157 > }
158  
159 <    if( have_cutoffgroups ){
160 <      sprintf( errMsg,
161 <               "MoleculeStamp error, n_cutoffgroups already declared for"
162 <               " molecule: %s\n",
177 <               name );
178 <      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      }
164 <    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;
164 >
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;
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 );
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;
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;
202 >    
203   }
204  
205 < char* MoleculeStamp::addAtom( AtomStamp* the_atom, int atomIndex ){
206 <  
207 <  if( have_atoms && atomIndex < n_atoms ) atoms[atomIndex] = the_atom;
208 <  else {
209 <    if( have_atoms ){
297 <      sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
298 <               atomIndex );
299 <      return strdup( errMsg );
205 > struct BendLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
206 >    bool operator()(IntTuple3 b1, IntTuple3 b2) {
207 >        return b1.first < b2.first
208 >             || (!(b2.first < b1.first) && b1.second < b2.second)
209 >             || (!(b2.first < b1.first) && !(b2.second < b2.second) && b1.third < b2.third);
210      }
211 <    else return strdup("MoleculeStamp error, nAtoms not given before"
302 <                       " first atom declaration." );
303 <  }
211 > };
212  
213 <  return NULL;
214 < }
215 <
216 < char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
217 <                                   int rigidBodyIndex ){
218 <  
219 <  if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
220 <    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 < }
213 > void MoleculeStamp::checkBends() {
214 >    for(int i = 0; i < getNBends(); ++i) {
215 >        BendStamp* bendStamp = getBendStamp(i);
216 >        std::vector<int> bendAtoms =  bendStamp->getMembers();
217 >        std::vector<int>::iterator j =std::find_if(bendAtoms.begin(), bendAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
218 >        if (j != bendAtoms.end()) {
219 >            std::cout << "Error in Molecule " << getName() << " : atoms of bend" << containerToString(bendAtoms) << "have invalid indices\n";
220 >        }
221  
222 < char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
223 <                                     int cutoffGroupIndex ){
224 <  
225 <  if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
226 <    cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
227 <  else {
228 <    if( have_cutoffgroups ){
229 <      sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
230 <               cutoffGroupIndex );
231 <      return strdup( errMsg );
222 >        if (bendAtoms.size() == 2 ) {
223 >            if (!bendStamp->haveGhostVectorSource()) {
224 >                std::cout << "Error in Molecule " << getName() << ": ghostVectorSouce is missing\n";
225 >            }else{
226 >                int ghostIndex = bendStamp->getGhostVectorSource();
227 >                if (ghostIndex < getNAtoms()) {
228 >                    if (std::find(bendAtoms.begin(), bendAtoms.end(), ghostIndex) == bendAtoms.end()) {
229 >                      std::cout <<  "Error in Molecule " << getName() << ": ghostVectorSouce "<< ghostIndex<<"is invalid\n";
230 >                    }
231 >                    if (!getAtomStamp(ghostIndex)->haveOrientation()) {
232 >                        std::cout <<  "Error in Molecule " << getName() << ": ghost atom must be a directioanl atom\n";
233 >                    }
234 >                }else {
235 >                    std::cout << "Error in Molecule " << getName() <<  ": ghostVectorsource " << ghostIndex<< "  is invalid\n";
236 >                }
237 >            }
238 >        } else if (bendAtoms.size() == 3 && bendStamp->haveGhostVectorSource()) {
239 >            std::cout <<  "Error in Molecule " << getName() << ": normal bend should not have ghostVectorSouce\n";
240 >        }
241      }
337 –    else return strdup("MoleculeStamp error, nCutoffGroups not given before"
338 –                       " first CutoffGroup declaration." );
339 –  }
340 –  
341 –  return NULL;
342 – }
242  
243 < char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
244 <  
245 <  
246 <  if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
247 <  else{
248 <    if( have_bonds ){
249 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
250 <               bondIndex );
251 <      return strdup( errMsg );
243 >    for(int i = 0; i < getNBends(); ++i) {
244 >        BendStamp* bendStamp = getBendStamp(i);
245 >        std::vector<int> bendAtoms =  bendStamp->getMembers();
246 >        std::vector<int> rigidSet(getNRigidBodies(), 0);
247 >        std::vector<int>::iterator j;
248 >        for( j = bendAtoms.begin(); j != bendAtoms.end(); ++j) {
249 >            int rigidbodyIndex = atom2Rigidbody[*j];
250 >            if (rigidbodyIndex >= 0) {
251 >                ++rigidSet[rigidbodyIndex];
252 >                if (rigidSet[rigidbodyIndex] > 1) {
253 >                    std::cout << "Error in Molecule " << getName() << ": bend" << containerToString(bendAtoms) << " belong to same rigidbody " << rigidbodyIndex << "\n";                    
254 >                }
255 >            }
256 >        }
257 >    }
258 >    
259 >    
260 >    std::set<IntTuple3, BendLessThan> allBends;
261 >    std::set<IntTuple3, BendLessThan>::iterator iter;
262 >    for(int i = 0; i < getNBends(); ++i) {
263 >        BendStamp* bendStamp= getBendStamp(i);
264 >        std::vector<int> bend = bendStamp->getMembers();
265 >        if (bend.size() == 2) {
266 >        // in case we have two ghost bend. For example,
267 >        // bend {
268 >        // members (0, 1);
269 >        //   ghostVectorSource = 0;
270 >        // }
271 >        // and
272 >        // bend {
273 >        //   members (0, 1);
274 >        // ghostVectorSource = 0;
275 >        // }
276 >        // In order to distinguish them. we expand them to Tuple3.
277 >        // the first one is expanded to (0, 0, 1) while the second one is expaned to (0, 1, 1)
278 >             int ghostIndex = bendStamp->getGhostVectorSource();
279 >             std::vector<int>::iterator j = std::find(bend.begin(), bend.end(), ghostIndex);
280 >             if (j != bend.end()) {
281 >                bend.insert(j, ghostIndex);
282 >             }
283 >        }
284 >        
285 >        IntTuple3 bendTuple(bend[0], bend[1], bend[2]);
286 >        //make sure bendTuple.first is always less than or equal to bendTuple.third
287 >        if (bendTuple.first > bendTuple.third) {
288 >            std::swap(bendTuple.first, bendTuple.third);
289 >        }
290 >        
291 >        iter = allBends.find(bendTuple);
292 >        if ( iter != allBends.end()) {
293 >            std::cout << "Error in Molecule " << getName() << ": " << "Bend appears multiple times\n";
294 >        } else {
295 >            allBends.insert(bendTuple);
296 >        }
297      }
354 –    else return strdup("MoleculeStamp error, nBonds not given before"
355 –                       "first bond declaration." );
356 –  }
298  
299 <  return NULL;
300 < }
299 >    for (int i = 0; i < getNBonds(); ++i) {
300 >        BondStamp* bondStamp = getBondStamp(i);
301 >        int a = bondStamp->getA();
302 >        int b = bondStamp->getB();
303  
304 < char* MoleculeStamp::addBend( BendStamp* the_bend, int bendIndex ){
305 <  
306 <  
307 <  if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
308 <  else{
309 <    if( have_bends ){
310 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
311 <               bendIndex );
312 <      return strdup( errMsg );
304 >        AtomStamp* atomA = getAtomStamp(a);
305 >        AtomStamp* atomB = getAtomStamp(b);
306 >
307 >        //find bend c--a--b
308 >        AtomStamp::AtomIter ai;
309 >        for(int c= atomA->getFirstBonedAtom(ai);c != -1;c = atomA->getNextBonedAtom(ai))
310 >        {
311 >            if(b == c)
312 >                continue;          
313 >            
314 >            IntTuple3 newBend(c, a, b);
315 >            if (newBend.first > newBend.third) {
316 >                std::swap(newBend.first, newBend.third);
317 >            }
318 >
319 >            if (allBends.find(newBend) == allBends.end() ) {                
320 >                allBends.insert(newBend);
321 >                BendStamp * newBendStamp = new BendStamp();
322 >                newBendStamp->setMembers(newBend);
323 >                addBendStamp(newBendStamp);
324 >            }
325 >        }        
326 >
327 >        //find bend a--b--c
328 >        for(int c= atomB->getFirstBonedAtom(ai);c != -1;c = atomB->getNextBonedAtom(ai))
329 >        {
330 >            if(a == c)
331 >                continue;          
332 >
333 >            IntTuple3 newBend( a, b, c);
334 >            if (newBend.first > newBend.third) {
335 >                std::swap(newBend.first, newBend.third);
336 >            }            
337 >            if (allBends.find(newBend) == allBends.end() ) {                
338 >                allBends.insert(newBend);
339 >                BendStamp * newBendStamp = new BendStamp();
340 >                newBendStamp->setMembers(newBend);
341 >                addBendStamp(newBendStamp);
342 >            }
343 >        }        
344      }
371 –    else return strdup("MoleculeStamp error, nBends not given before"
372 –                       "first bend declaration." );
373 –  }
345  
375 –  return NULL;
346   }
347  
348 < char* MoleculeStamp::addTorsion( TorsionStamp* the_torsion, int torsionIndex ){
349 <  
350 <  
351 <  if( have_torsions && torsionIndex < n_torsions )
352 <    torsions[torsionIndex] = the_torsion;
353 <  else{
354 <    if( have_torsions ){
385 <      sprintf( errMsg, "MoleculeStamp error, %d out of nTorsions range",
386 <               torsionIndex );
387 <      return strdup( errMsg );
348 > struct TorsionLessThan : public std::binary_function<IntTuple4, IntTuple4, bool> {
349 >    bool operator()(IntTuple4 t1, IntTuple4 t2) {
350 >
351 >        return t1.first < t2.first
352 >             || (!(t2.first < t1.first) && t1.second < t2.second)
353 >             || (!(t2.first < t1.first) && !(t2.second < t2.second) && t1.third < t2.third)
354 >             ||(!(t2.first < t1.first) && !(t2.second < t2.second) && !(t2.third < t1.third) && t1.fourth < t2.fourth);
355      }
389 –    else return strdup("MoleculeStamp error, nTorsions not given before"
390 –                       "first torsion declaration." );
391 –  }
356  
393 –  return NULL;
394 – }
357  
396 – char* MoleculeStamp::checkMe( void ){
397 –  
398 –  int i;
399 –  short int no_atom, no_rigidbody, no_cutoffgroup;
358  
359 <  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 <  }
359 > };
360  
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 –  }
361  
362 <  no_cutoffgroup = 0;
363 <  for( i=0; i<n_cutoffgroups; i++ ){
364 <    if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
365 <  }
362 > void MoleculeStamp::checkTorsions() {
363 >    for(int i = 0; i < getNTorsions(); ++i) {
364 >        TorsionStamp* torsionStamp = getTorsionStamp(i);
365 >        std::vector<int> torsionAtoms =  torsionStamp ->getMembers();
366 >        std::vector<int>::iterator j =std::find_if(torsionAtoms.begin(), torsionAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
367 >        if (j != torsionAtoms.end()) {
368 >            std::cout << "Error in Molecule " << getName() << ": atoms of torsion" << containerToString(torsionAtoms) << " have invalid indices\n";
369 >        }
370 >    }
371 >    
372 >    for(int i = 0; i < getNTorsions(); ++i) {
373 >        TorsionStamp* torsionStamp = getTorsionStamp(i);
374 >        std::vector<int> torsionAtoms =  torsionStamp->getMembers();
375 >        std::vector<int> rigidSet(getNRigidBodies(), 0);
376 >        std::vector<int>::iterator j;
377 >        for( j = torsionAtoms.begin(); j != torsionAtoms.end(); ++j) {
378 >            int rigidbodyIndex = atom2Rigidbody[*j];
379 >            if (rigidbodyIndex >= 0) {
380 >                ++rigidSet[rigidbodyIndex];
381 >                if (rigidSet[rigidbodyIndex] > 1) {
382 >                    std::cout << "Error in Molecule " << getName() << ": torsion" << containerToString(torsionAtoms) << "is invalid\n";                  
383 >                }
384 >            }
385 >        }
386 >    }    
387  
388 <  if( no_cutoffgroup ){
389 <    sprintf( errMsg,
390 <             "MoleculeStamp error. Not all of the CutoffGroups were"
391 <             " declared in molecule \"%s\".\n", name );
392 <    return strdup( errMsg );
393 <  }
394 <  
395 <  no_atom = 0;
396 <  for( i=0; i<n_atoms; i++ ){
397 <    if( atoms[i] == NULL ) no_atom = 1;
398 <  }
388 >    std::set<IntTuple4, TorsionLessThan> allTorsions;
389 >    std::set<IntTuple4, TorsionLessThan>::iterator iter;
390 >     for(int i = 0; i < getNTorsions(); ++i) {
391 >         TorsionStamp* torsionStamp= getTorsionStamp(i);
392 >         std::vector<int> torsion = torsionStamp->getMembers();
393 >         if (torsion.size() == 3) {
394 >             int ghostIndex = torsionStamp->getGhostVectorSource();
395 >             std::vector<int>::iterator j = std::find(torsion.begin(), torsion.end(), ghostIndex);
396 >             if (j != torsion.end()) {
397 >                torsion.insert(j, ghostIndex);
398 >             }
399 >         }
400  
401 <  if( no_atom ){
402 <    sprintf( errMsg,
403 <             "MoleculeStamp error. Not all of the atoms were"
404 <             " declared in molecule \"%s\".\n", name );
405 <    return strdup( errMsg );
442 <  }
401 >        IntTuple4 torsionTuple(torsion[0], torsion[1], torsion[2], torsion[3]);
402 >        if (torsionTuple.first > torsionTuple.fourth) {
403 >            std::swap(torsionTuple.first, torsionTuple.fourth);
404 >            std::swap(torsionTuple.second, torsionTuple.third);                    
405 >        }                
406  
407 <  n_integrable = n_atoms;
408 <  for (i = 0; i < n_rigidbodies; i++)
409 <    n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
410 <  
411 <  if (n_integrable <= 0 || n_integrable > n_atoms) {
412 <    sprintf( errMsg,
413 <             "MoleculeStamp error. n_integrable is either <= 0 or"
451 <             " greater than n_atoms in molecule \"%s\".\n", name );
452 <    return strdup( errMsg );
453 <  }
407 >         iter = allTorsions.find(torsionTuple);
408 >         if ( iter == allTorsions.end()) {
409 >            allTorsions.insert(torsionTuple);
410 >         } else {
411 >            std::cout << "Error in Molecule " << getName() << ": " << "Torsion appears multiple times\n";
412 >         }
413 >     }
414  
415 <  return NULL;
416 < }  
415 >    for (int i = 0; i < getNBonds(); ++i) {
416 >        BondStamp* bondStamp = getBondStamp(i);
417 >        int b = bondStamp->getA();
418 >        int c = bondStamp->getB();
419  
420 +        AtomStamp* atomB = getAtomStamp(b);
421 +        AtomStamp* atomC = getAtomStamp(c);
422  
423 +        AtomStamp::AtomIter ai2;
424 +        AtomStamp::AtomIter ai3;
425 +
426 +        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
427 +        {
428 +            if(a == c)
429 +                continue;
430 +
431 +            for(int d = atomC->getFirstBonedAtom(ai3);d != -1;d = atomC->getNextBonedAtom(ai3))
432 +            {
433 +                if(d == b)
434 +                    continue;
435 +                
436 +                IntTuple4 newTorsion(a, b, c, d);
437 +                //make sure the first element is always less than or equal to the fourth element in IntTuple4
438 +                if (newTorsion.first > newTorsion.fourth) {
439 +                    std::swap(newTorsion.first, newTorsion.fourth);
440 +                    std::swap(newTorsion.second, newTorsion.third);                    
441 +                }                
442 +                if (allTorsions.find(newTorsion) == allTorsions.end() ) {                
443 +                    allTorsions.insert(newTorsion);
444 +                    TorsionStamp * newTorsionStamp = new TorsionStamp();
445 +                    newTorsionStamp->setMembers(newTorsion);
446 +                    addTorsionStamp(newTorsionStamp);                    
447 +                }            
448 +            }
449 +        }    
450 +    }
451 +
452 + }
453 +
454 + void MoleculeStamp::checkRigidBodies() {
455 +     std::vector<RigidBodyStamp*>::iterator ri = std::find(rigidBodyStamps_.begin(), rigidBodyStamps_.end(), static_cast<RigidBodyStamp*>(NULL));
456 +     if (ri != rigidBodyStamps_.end()) {
457 +         std::cout << "Error in Molecule " << getName() << ":rigidBody[" <<  ri - rigidBodyStamps_.begin()<< "] is missing\n";
458 +     }
459 +
460 +    for (int i = 0; i < getNRigidBodies(); ++i) {
461 +        RigidBodyStamp* rbStamp = getRigidBodyStamp(i);
462 +        std::vector<int> rigidAtoms =  rbStamp ->getMembers();
463 +        std::vector<int>::iterator j =std::find_if(rigidAtoms.begin(), rigidAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
464 +        if (j != rigidAtoms.end()) {
465 +            std::cout << "Error in Molecule " << getName();
466 +        }
467 +        
468 +    }    
469 + }
470 +
471 + void MoleculeStamp::checkCutoffGroups() {
472 +
473 +    for(int i = 0; i < getNCutoffGroups(); ++i) {
474 +        CutoffGroupStamp* cutoffGroupStamp = getCutoffGroupStamp(i);
475 +        std::vector<int> cutoffGroupAtoms =  cutoffGroupStamp ->getMembers();
476 +        std::vector<int>::iterator j =std::find_if(cutoffGroupAtoms.begin(), cutoffGroupAtoms.end(), std::bind2nd(std::greater<int>(), getNAtoms()-1));
477 +        if (j != cutoffGroupAtoms.end()) {
478 +            std::cout << "Error in Molecule " << getName() << ": cutoffGroup" << " is out of range\n";
479 +        }
480 +    }    
481 + }
482 +
483 + void MoleculeStamp::checkFragments() {
484 +
485 +    std::vector<FragmentStamp*>::iterator fi = std::find(fragmentStamps_.begin(), fragmentStamps_.end(), static_cast<FragmentStamp*>(NULL));
486 +    if (fi != fragmentStamps_.end()) {
487 +        std::cout << "Error in Molecule " << getName() << ":fragment[" <<  fi - fragmentStamps_.begin()<< "] is missing\n";
488 +    }
489 +    
490 + }
491 +
492 + void MoleculeStamp::fillBondInfo() {
493 +
494 +    for (int i = 0; i < getNBonds(); ++i) {
495 +        BondStamp* bondStamp = getBondStamp(i);
496 +        int a = bondStamp->getA();
497 +        int b = bondStamp->getB();
498 +        AtomStamp* atomA = getAtomStamp(a);
499 +        AtomStamp* atomB = getAtomStamp(b);
500 +        atomA->addBond(i);
501 +        atomA->addBondedAtom(b);
502 +        atomB->addBond(i);        
503 +        atomB->addBondedAtom(a);
504 +
505 +    }
506 + }
507 +
508 +
509 +
510   //Function Name: isBondInSameRigidBody
511   //Return true is both atoms of the bond belong to the same rigid body, otherwise return false
512   bool MoleculeStamp::isBondInSameRigidBody(BondStamp* bond){
# Line 479 | Line 530 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
530   // Function Name: isAtomInRigidBody
531   //return false if atom does not belong to a rigid body, otherwise return true
532   bool MoleculeStamp::isAtomInRigidBody(int atomIndex){
533 <  int whichRigidBody;
483 <  int consAtomIndex;
484 <
485 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
533 >  return atom2Rigidbody[atomIndex] >=0 ;
534    
535   }
536  
# Line 493 | Line 541 | bool MoleculeStamp::isAtomInRigidBody(int atomIndex, i
541   //whichRigidBody: the index of rigidbody in component
542   //consAtomIndex:  the position of joint atom apears in  rigidbody's definition
543   bool MoleculeStamp::isAtomInRigidBody(int atomIndex, int& whichRigidBody, int& consAtomIndex){
496 –  RigidBodyStamp* rbStamp;
497 –  int numRb;
498 –  int numAtom;
544  
545 +  
546 +
547    whichRigidBody = -1;
548    consAtomIndex = -1;
549  
550 <  numRb = this->getNRigidBodies();
551 <  
552 <  for(int i = 0 ; i < numRb; i++){
553 <    rbStamp = this->getRigidBody(i);
554 <    numAtom = rbStamp->getNMembers();
555 <    for(int j = 0; j < numAtom; j++)
509 <      if (rbStamp->getMember(j) == atomIndex){
510 <        whichRigidBody = i;
550 >  if (atom2Rigidbody[atomIndex] >=0) {
551 >    whichRigidBody = atom2Rigidbody[atomIndex];
552 >    RigidBodyStamp* rbStamp = getRigidBodyStamp(whichRigidBody);
553 >    int numAtom = rbStamp->getNMembers();
554 >    for(int j = 0; j < numAtom; j++) {
555 >      if (rbStamp->getMemberAt(j) == atomIndex){
556          consAtomIndex = j;
557          return true;
558        }
559 +    }
560    }
561  
562    return false;
# Line 520 | Line 566 | vector<pair<int, int> > MoleculeStamp::getJointAtoms(i
566   //return the position of joint atom apears in  rigidbody's definition
567   //for the time being, we will use the most inefficient algorithm, the complexity is O(N2)
568   //actually we could improve the complexity to O(NlgN) by sorting the atom index in rigid body first
569 < vector<pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
569 > std::vector<std::pair<int, int> > MoleculeStamp::getJointAtoms(int rb1, int rb2){
570    RigidBodyStamp* rbStamp1;
571    RigidBodyStamp* rbStamp2;
572    int natomInRb1;
573    int natomInRb2;
574    int atomIndex1;
575    int atomIndex2;
576 <  vector<pair<int, int> > jointAtomIndexPair;
576 >  std::vector<std::pair<int, int> > jointAtomIndexPair;
577    
578 <  rbStamp1 = this->getRigidBody(rb1);
578 >  rbStamp1 = this->getRigidBodyStamp(rb1);
579    natomInRb1 =rbStamp1->getNMembers();
580  
581 <  rbStamp2 = this->getRigidBody(rb2);
581 >  rbStamp2 = this->getRigidBodyStamp(rb2);
582    natomInRb2 =rbStamp2->getNMembers();
583  
584    for(int i = 0; i < natomInRb1; i++){
585 <    atomIndex1 = rbStamp1->getMember(i);
585 >    atomIndex1 = rbStamp1->getMemberAt(i);
586        
587      for(int j= 0; j < natomInRb1; j++){
588 <      atomIndex2 = rbStamp2->getMember(j);
588 >      atomIndex2 = rbStamp2->getMemberAt(j);
589  
590        if(atomIndex1 == atomIndex2){
591 <        jointAtomIndexPair.push_back(make_pair(i, j));
591 >        jointAtomIndexPair.push_back(std::make_pair(i, j));
592          break;
593        }
594        
595 <    }//end for(j =0)
595 >    }
596  
597 <  }//end for (i = 0)
597 >  }
598  
599    return jointAtomIndexPair;
600   }
601 +
602 + }

Diff Legend

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