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 2507 by tim, Sat Dec 10 16:54:40 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" ) ){
139 <    n_bends = (int)rhs;
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 >    }
140  
141 <    if( have_bends ){
142 <      sprintf( errMsg,
143 <               "MoleculeStamp error, n_bends already declared for"
144 <               " molecule: %s\n",
145 <               name);
146 <      return strdup( errMsg );
147 <    }
148 <    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;
141 >    checkAtoms();
142 >    checkBonds();
143 >    fillBondInfo();
144 >    checkBends();
145 >    checkTorsions();
146 >    checkRigidBodies();
147 >    checkCutoffGroups();
148 >    checkFragments();
149  
150 <    if( have_torsions ){
151 <      sprintf( errMsg,
152 <               "MoleculeStamp error, n_torsions already declared for"
153 <               " molecule: %s\n",
147 <               name );
148 <      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_torsions = 1;
151 <    torsions = new TorsionStamp*[n_torsions];
152 <    for( i=0; i<n_torsions; i++ ) torsions[i] = NULL;
153 <  }
155 >    nintegrable_ = getNAtoms()+ getNRigidBodies() - nrigidAtoms;
156  
157 <  else if( !strcmp( lhs, "nRigidBodies" ) ){
156 <    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",
162 <               name );
163 <      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      }
165 –    have_rigidbodies = 1;
166 –    rigidBodies = new RigidBodyStamp*[n_rigidbodies];
167 –    for( i=0; i<n_rigidbodies; i++ ) rigidBodies[i] = NULL;
168 –  }
164  
170 –  else if( !strcmp( lhs, "nCutoffGroups" ) ){
171 –    n_cutoffgroups = (int)rhs;
172 –
173 –    if( have_cutoffgroups ){
174 –      sprintf( errMsg,
175 –               "MoleculeStamp error, n_cutoffgroups already declared for"
176 –               " molecule: %s\n",
177 –               name );
178 –      return strdup( errMsg );
179 –    }
180 –    have_cutoffgroups = 1;
181 –    cutoffGroups = new CutoffGroupStamp*[n_cutoffgroups];
182 –    for( i=0; i<n_cutoffgroups; i++ ) cutoffGroups[i] = NULL;
183 –  }
184 –  
185 –  else{
186 –    if( unhandled == NULL ) unhandled = new LinkedAssign( lhs, rhs );
187 –    else unhandled->add( lhs, rhs );
188 –    have_extras = 1;
189 –  }
190 –  return NULL;
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 ){
210 <      sprintf( errMsg, "MoleculeStamp error, %d out of nAtoms range",
211 <               atomIndex );
212 <      return strdup( errMsg );
300 <    }
301 <    else return strdup("MoleculeStamp error, nAtoms not given before"
302 <                       " first atom declaration." );
303 <  }
304 <
305 <  return NULL;
306 < }
307 <
308 < char* MoleculeStamp::addRigidBody( RigidBodyStamp* the_rigidbody,
309 <                                   int rigidBodyIndex ){
310 <  
311 <  if( have_rigidbodies && rigidBodyIndex < n_rigidbodies )
312 <    rigidBodies[rigidBodyIndex] = the_rigidbody;
313 <  else {
314 <    if( have_rigidbodies ){
315 <      sprintf( errMsg, "MoleculeStamp error, %d out of nRigidBodies range",
316 <               rigidBodyIndex );
317 <      return strdup( errMsg );
318 <    }
319 <    else return strdup("MoleculeStamp error, nRigidBodies not given before"
320 <                       " first rigidBody declaration." );
321 <  }
322 <  
323 <  return NULL;
324 < }
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 < char* MoleculeStamp::addCutoffGroup( CutoffGroupStamp* the_cutoffgroup,
215 <                                     int cutoffGroupIndex ){
216 <  
217 <  if( have_cutoffgroups && cutoffGroupIndex < n_cutoffgroups )
218 <    cutoffGroups[cutoffGroupIndex] = the_cutoffgroup;
219 <  else {
220 <    if( have_cutoffgroups ){
221 <      sprintf( errMsg, "MoleculeStamp error, %d out of nCutoffGroups range",
222 <               cutoffGroupIndex );
223 <      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      }
337 –    else return strdup("MoleculeStamp error, nCutoffGroups not given before"
338 –                       " first CutoffGroup declaration." );
339 –  }
340 –  
341 –  return NULL;
342 – }
234  
235 < char* MoleculeStamp::addBond( BondStamp* the_bond, int bondIndex ){
236 <  
237 <  
238 <  if( have_bonds && bondIndex < n_bonds ) bonds[bondIndex] = the_bond;
239 <  else{
240 <    if( have_bonds ){
241 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBonds range",
242 <               bondIndex );
243 <      return strdup( errMsg );
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      }
354 –    else return strdup("MoleculeStamp error, nBonds not given before"
355 –                       "first bond declaration." );
356 –  }
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::addBend( BendStamp* the_bend, int bendIndex ){
297 <  
298 <  
299 <  if( have_bends && bendIndex < n_bends ) bends[bendIndex] = the_bend;
300 <  else{
301 <    if( have_bends ){
302 <      sprintf( errMsg, "MoleculeStamp error, %d out of nBends range",
303 <               bendIndex );
304 <      return strdup( errMsg );
296 >        AtomStamp* atomA = getAtomStamp(a);
297 >        AtomStamp* atomB = getAtomStamp(b);
298 >
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 >            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 >        //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 >            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      }
371 –    else return strdup("MoleculeStamp error, nBends not given before"
372 –                       "first bend declaration." );
373 –  }
337  
375 –  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",
386 <               torsionIndex );
387 <      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;
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 >        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 >         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 >    for (int i = 0; i < getNBonds(); ++i) {
394 >        BondStamp* bondStamp = getBondStamp(i);
395 >        int b = bondStamp->getA();
396 >        int c = bondStamp->getB();
397 >
398 >        AtomStamp* atomB = getAtomStamp(b);
399 >        AtomStamp* atomC = getAtomStamp(c);
400 >
401 >        AtomStamp::AtomIter ai2;
402 >        AtomStamp::AtomIter ai3;
403 >
404 >        for(int a = atomB->getFirstBonedAtom(ai2);a != -1;a = atomB->getNextBonedAtom(ai2))
405 >        {
406 >            if(a == c)
407 >                continue;
408 >
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   }
431  
432 < char* MoleculeStamp::checkMe( void ){
433 <  
434 <  int i;
435 <  short int no_atom, no_rigidbody, no_cutoffgroup;
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 <  if( !have_name ) return strdup( "MoleculeStamp error. Molecule's name"
439 <                                  " was not given.\n" );
440 <  
441 <  if( !have_atoms ){
442 <    return strdup( "MoleculeStamp error. Molecule contains no atoms." );
443 <  }
444 <  
445 <  no_rigidbody = 0;
446 <  for( i=0; i<n_rigidbodies; i++ ){
447 <    if( rigidBodies[i] == NULL ) no_rigidbody = 1;
411 <  }
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 <  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 <  }
449 > void MoleculeStamp::checkCutoffGroups() {
450  
451 <  no_cutoffgroup = 0;
452 <  for( i=0; i<n_cutoffgroups; i++ ){
453 <    if( cutoffGroups[i] == NULL ) no_cutoffgroup = 1;
454 <  }
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 <  if( no_cutoffgroup ){
426 <    sprintf( errMsg,
427 <             "MoleculeStamp error. Not all of the CutoffGroups were"
428 <             " declared in molecule \"%s\".\n", name );
429 <    return strdup( errMsg );
430 <  }
431 <  
432 <  no_atom = 0;
433 <  for( i=0; i<n_atoms; i++ ){
434 <    if( atoms[i] == NULL ) no_atom = 1;
435 <  }
461 > void MoleculeStamp::checkFragments() {
462  
463 <  if( no_atom ){
464 <    sprintf( errMsg,
465 <             "MoleculeStamp error. Not all of the atoms were"
466 <             " declared in molecule \"%s\".\n", name );
467 <    return strdup( errMsg );
468 <  }
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 <  n_integrable = n_atoms;
445 <  for (i = 0; i < n_rigidbodies; i++)
446 <    n_integrable = n_integrable - rigidBodies[i]->getNMembers() + 1; //rigidbody is an integrable object
447 <  
448 <  if (n_integrable <= 0 || n_integrable > n_atoms) {
449 <    sprintf( errMsg,
450 <             "MoleculeStamp error. n_integrable is either <= 0 or"
451 <             " greater than n_atoms in molecule \"%s\".\n", name );
452 <    return strdup( errMsg );
453 <  }
470 > void MoleculeStamp::fillBondInfo() {
471  
472 <  return NULL;
473 < }  
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 479 | 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;
483 <  int consAtomIndex;
484 <
485 <  return isAtomInRigidBody(atomIndex, whichRigidBody, consAtomIndex);
511 >  return atom2Rigidbody[atomIndex] >=0 ;
512    
513   }
514  
# Line 493 | 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){
496 –  RigidBodyStamp* rbStamp;
497 –  int numRb;
498 –  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++)
509 <      if (rbStamp->getMember(j) == atomIndex){
510 <        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 520 | Line 544 | vector<pair<int, int> > MoleculeStamp::getJointAtoms(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 < vector<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 <  vector<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(make_pair(i, j));
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