ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp
(Generate patch)

Comparing:
trunk/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1735 by tim, Fri Nov 12 17:40:03 2004 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
2 < #include <string.h>
3 < #include <math.h>
1 > /*
2 > * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 > *
4 > * Contact: oopse@oopse.org
5 > *
6 > * This program is free software; you can redistribute it and/or
7 > * modify it under the terms of the GNU Lesser General Public License
8 > * as published by the Free Software Foundation; either version 2.1
9 > * of the License, or (at your option) any later version.
10 > * All we ask is that proper credit is given for our work, which includes
11 > * - but is not limited to - adding the above copyright notice to the beginning
12 > * of your source code files, and to any copyright notice that you may distribute
13 > * with programs based on this work.
14 > *
15 > * This program is distributed in the hope that it will be useful,
16 > * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 > * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 > * GNU Lesser General Public License for more details.
19 > *
20 > * You should have received a copy of the GNU Lesser General Public License
21 > * along with this program; if not, write to the Free Software
22 > * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23 > *
24 > */
25  
26 < #include <iostream>
27 < using namespace std;
26 > /**
27 > * @file SimInfo.cpp
28 > * @author    tlin
29 > * @date  11/02/2004
30 > * @version 1.0
31 > */
32  
33 < #include "SimInfo.hpp"
9 < #define __C
10 < #include "fSimulation.h"
11 < #include "simError.h"
33 > #include <algorithm>
34  
35 < #include "fortranWrappers.hpp"
35 > #include "brains/SimInfo.hpp"
36 > #include "utils/MemoryUtils.hpp"
37  
38 < #include "MatVec3.h"
38 > namespace oopse {
39  
40 < #ifdef IS_MPI
41 < #include "mpiSimulation.hpp"
42 < #endif
40 > SimInfo::SimInfo(const std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
41 >                                ForceField* ff, Globals* globals) :
42 >                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
43 >                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
44 >                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) {
45  
46 < inline double roundMe( double x ){
47 <  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
48 < }
49 <          
50 < inline double min( double a, double b ){
51 <  return (a < b ) ? a : b;
52 < }
46 >    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
47 >    int nCutoffAtoms; // number of atoms belong to cutoff groups
48 >    int ngroups;          //total cutoff groups defined in meta-data file
49 >    MoleculeStamp* molStamp;
50 >    int nMolWithSameStamp;
51 >    CutoffGroupStamp* cgStamp;
52 >    int nAtomsIngroups;
53 >    int nCutoffGroupsInStamp;    
54  
55 < SimInfo* currentInfo;
55 >    nGlobalAtoms_ =  0;
56 >    ngroups = 0;
57 >    
58 >    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
59 >        molStamp = i->first;
60 >        nMolWithSameStamp = i->second;
61 >        
62 >        addMoleculeStamp(molStamp, nMolWithSameStamp);
63 >        
64 >        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
65 >        
66 >        nAtomsIngroups = 0;
67 >        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
68 >        
69 >        for (int j=0; j < nCutoffGroupsInStamp; j++) {
70 >            cgStamp = molStamp->getCutoffGroup(j);
71 >            nAtomsIngroups += cgStamp->getNMembers();
72 >        }
73  
74 < SimInfo::SimInfo(){
74 >        ngroups += *nMolWithSameStamp;
75 >        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
76 >    }
77  
78 <  n_constraints = 0;
79 <  nZconstraints = 0;
80 <  n_oriented = 0;
81 <  n_dipoles = 0;
82 <  ndf = 0;
38 <  ndfRaw = 0;
39 <  nZconstraints = 0;
40 <  the_integrator = NULL;
41 <  setTemp = 0;
42 <  thermalTime = 0.0;
43 <  currentTime = 0.0;
44 <  rCut = 0.0;
45 <  rSw = 0.0;
78 >    //every free atom (atom does not belong to cutoff groups) is a cutoff group
79 >    //therefore the total number of cutoff groups in the system is equal to
80 >    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
81 >    //file plus the number of cutoff groups defined in meta-data file
82 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups;
83  
84 <  haveRcut = 0;
85 <  haveRsw = 0;
49 <  boxIsInit = 0;
50 <  
51 <  resetTime = 1e99;
84 >    //initialize globalGroupMembership_, every element of this array will be 0
85 >    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
86  
87 <  orthoRhombic = 0;
54 <  orthoTolerance = 1E-6;
55 <  useInitXSstate = true;
87 >    nGlobalMols_ = molStampIds_.size();
88  
89 <  usePBC = 0;
90 <  useLJ = 0;
91 <  useSticky = 0;
92 <  useCharges = 0;
93 <  useDipoles = 0;
62 <  useReactionField = 0;
63 <  useGB = 0;
64 <  useEAM = 0;
65 <  useSolidThermInt = 0;
66 <  useLiquidThermInt = 0;
89 > #ifdef IS_MPI    
90 >    molToProcMap_.resize(nGlobalMols_);
91 > #endif
92 >    
93 > }
94  
95 <  haveCutoffGroups = false;
95 > SimInfo::~SimInfo() {
96 >    //MemoryUtils::deleteVectorOfPointer(molecules_);
97  
98 <  excludes = Exclude::Instance();
98 >    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
99 >    
100 >    delete sman_;
101 >    delete globals_;
102 >    delete forceField_;
103  
72  myConfiguration = new SimState();
73
74  has_minimizer = false;
75  the_minimizer =NULL;
76
77  ngroup = 0;
78
79  wrapMeSimInfo( this );
104   }
105  
106  
107 < SimInfo::~SimInfo(){
107 > bool SimInfo::addMolecule(Molecule* mol) {
108 >    MoleculeIterator i;
109  
110 <  delete myConfiguration;
110 >    i = molecules_.find(mol->getGlobalIndex());
111 >    if (i != molecules_.end() ) {
112  
113 <  map<string, GenericData*>::iterator i;
114 <  
115 <  for(i = properties.begin(); i != properties.end(); i++)
116 <    delete (*i).second;
113 >        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
114 >        
115 >        nAtoms_ += mol->getNAtoms();
116 >        nBonds_ += mol->getNBonds();
117 >        nBends_ += mol->getNBends();
118 >        nTorsions_ += mol->getNTorsions();
119 >        nRigidBodies_ += mol->getNRigidBodies();
120 >        nIntegrableObjects_ += mol->getNIntegrableObjects();
121 >        nCutoffGroups_ += mol->getNCutoffGroups();
122 >        nConstraints_ += mol->getNConstraints();
123  
124 +        return true;
125 +    } else {
126 +        return false;
127 +    }
128   }
129  
130 < void SimInfo::setBox(double newBox[3]) {
131 <  
132 <  int i, j;
97 <  double tempMat[3][3];
130 > bool SimInfo::removeMolecule(Molecule* mol) {
131 >    MoleculeIterator i;
132 >    i = molecules_.find(mol->getGlobalIndex());
133  
134 <  for(i=0; i<3; i++)
100 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
134 >    if (i != molecules_.end() ) {
135  
136 <  tempMat[0][0] = newBox[0];
137 <  tempMat[1][1] = newBox[1];
138 <  tempMat[2][2] = newBox[2];
136 >        assert(mol == i->second);
137 >        
138 >        nAtoms_ -= mol->getNAtoms();
139 >        nBonds_ -= mol->getNBonds();
140 >        nBends_ -= mol->getNBends();
141 >        nTorsions_ -= mol->getNTorsions();
142 >        nRigidBodies_ -= mol->getNRigidBodies();
143 >        nIntegrableObjects_ -= mol->getNIntegrableObjects();
144 >        nCutoffGroups_ -= mol->getNCutoffGroups();
145 >        nConstraints_ -= mol->getNConstraints();
146  
147 <  setBoxM( tempMat );
147 >        molecules_.erase(mol->getGlobalIndex());
148  
149 < }
149 >        delete mol;
150 >        
151 >        return true;
152 >    } else {
153 >        return false;
154 >    }
155  
110 void SimInfo::setBoxM( double theBox[3][3] ){
111  
112  int i, j;
113  double FortranHmat[9]; // to preserve compatibility with Fortran the
114                         // ordering in the array is as follows:
115                         // [ 0 3 6 ]
116                         // [ 1 4 7 ]
117                         // [ 2 5 8 ]
118  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
156  
157 <  if( !boxIsInit ) boxIsInit = 1;
157 > }    
158  
159 <  for(i=0; i < 3; i++)
160 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
161 <  
162 <  calcBoxL();
163 <  calcHmatInv();
159 >        
160 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
161 >    i = molecules_.begin();
162 >    return i == molecules_.end() ? NULL : *i;
163 > }    
164  
165 <  for(i=0; i < 3; i++) {
166 <    for (j=0; j < 3; j++) {
167 <      FortranHmat[3*j + i] = Hmat[i][j];
131 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
132 <    }
133 <  }
134 <
135 <  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
136 <
165 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
166 >    ++i;
167 >    return i == molecules_.end() ? NULL : *i;    
168   }
138
169  
140 void SimInfo::getBoxM (double theBox[3][3]) {
170  
171 <  int i, j;
172 <  for(i=0; i<3; i++)
173 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
174 < }
171 > void SimInfo::calcNdf() {
172 >    int ndf_local;
173 >    MoleculeIterator i;
174 >    std::vector<StuntDouble*>::iterator j;
175 >    Molecule* mol;
176 >    StuntDouble* integrableObject;
177  
178 +    ndf_local = 0;
179 +    
180 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
181 +        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
182 +               integrableObject = mol->nextIntegrableObject(j)) {
183  
184 < void SimInfo::scaleBox(double scale) {
149 <  double theBox[3][3];
150 <  int i, j;
184 >            ndf_local += 3;
185  
186 <  // cerr << "Scaling box by " << scale << "\n";
186 >            if (integrableObject->isDirectional()) {
187 >                if (integrableObject->isLinear()) {
188 >                    ndf_local += 2;
189 >                } else {
190 >                    ndf_local += 3;
191 >                }
192 >            }
193 >            
194 >        }//end for (integrableObject)
195 >    }// end for (mol)
196 >    
197 >    // n_constraints is local, so subtract them on each processor
198 >    ndf_local -= nConstraints_;
199  
200 <  for(i=0; i<3; i++)
201 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
200 > #ifdef IS_MPI
201 >    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
202 > #else
203 >    ndf_ = ndf_local;
204 > #endif
205  
206 <  setBoxM(theBox);
206 >    // nZconstraints_ is global, as are the 3 COM translations for the
207 >    // entire system:
208 >    ndf_ = ndf_ - 3 - nZconstraints_;
209  
210   }
211  
212 < void SimInfo::calcHmatInv( void ) {
213 <  
163 <  int oldOrtho;
164 <  int i,j;
165 <  double smallDiag;
166 <  double tol;
167 <  double sanity[3][3];
212 > void SimInfo::calcNdfRaw() {
213 >    int ndfRaw_local;
214  
215 <  invertMat3( Hmat, HmatInv );
215 >    MoleculeIterator i;
216 >    std::vector<StuntDouble*>::iterator j;
217 >    Molecule* mol;
218 >    StuntDouble* integrableObject;
219  
220 <  // check to see if Hmat is orthorhombic
221 <  
222 <  oldOrtho = orthoRhombic;
223 <
224 <  smallDiag = fabs(Hmat[0][0]);
225 <  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
177 <  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
178 <  tol = smallDiag * orthoTolerance;
220 >    // Raw degrees of freedom that we have to set
221 >    ndfRaw_local = 0;
222 >    
223 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
224 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
225 >               integrableObject = mol->nextIntegrableObject(j)) {
226  
227 <  orthoRhombic = 1;
181 <  
182 <  for (i = 0; i < 3; i++ ) {
183 <    for (j = 0 ; j < 3; j++) {
184 <      if (i != j) {
185 <        if (orthoRhombic) {
186 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
187 <        }        
188 <      }
189 <    }
190 <  }
227 >            ndfRaw_local += 3;
228  
229 <  if( oldOrtho != orthoRhombic ){
230 <    
231 <    if( orthoRhombic ) {
232 <      sprintf( painCave.errMsg,
233 <               "OOPSE is switching from the default Non-Orthorhombic\n"
234 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
235 <               "\tThis is usually a good thing, but if you wan't the\n"
236 <               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
237 <               "\tvariable ( currently set to %G ) smaller.\n",
201 <               orthoTolerance);
202 <      painCave.severity = OOPSE_INFO;
203 <      simError();
229 >            if (integrableObject->isDirectional()) {
230 >                if (integrableObject->isLinear()) {
231 >                    ndfRaw_local += 2;
232 >                } else {
233 >                    ndfRaw_local += 3;
234 >                }
235 >            }
236 >            
237 >        }
238      }
239 <    else {
240 <      sprintf( painCave.errMsg,
241 <               "OOPSE is switching from the faster Orthorhombic to the more\n"
242 <               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
243 <               "\tThis is usually because the box has deformed under\n"
244 <               "\tNPTf integration. If you wan't to live on the edge with\n"
211 <               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
212 <               "\tvariable ( currently set to %G ) larger.\n",
213 <               orthoTolerance);
214 <      painCave.severity = OOPSE_WARNING;
215 <      simError();
216 <    }
217 <  }
239 >    
240 > #ifdef IS_MPI
241 >    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
242 > #else
243 >    ndfRaw_ = ndfRaw_local;
244 > #endif
245   }
246  
247 < void SimInfo::calcBoxL( void ){
247 > void SimInfo::calcNdfTrans() {
248 >    int ndfTrans_local;
249  
250 <  double dx, dy, dz, dsq;
250 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
251  
224  // boxVol = Determinant of Hmat
252  
253 <  boxVol = matDet3( Hmat );
253 > #ifdef IS_MPI
254 >    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
255 > #else
256 >    ndfTrans_ = ndfTrans_local;
257 > #endif
258  
259 <  // boxLx
260 <  
261 <  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
231 <  dsq = dx*dx + dy*dy + dz*dz;
232 <  boxL[0] = sqrt( dsq );
233 <  //maxCutoff = 0.5 * boxL[0];
259 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraints_;
260 >
261 > }
262  
263 <  // boxLy
264 <  
265 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
266 <  dsq = dx*dx + dy*dy + dz*dz;
267 <  boxL[1] = sqrt( dsq );
268 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
263 > void SimInfo::addExcludePairs(Molecule* mol) {
264 >    std::vector<Bond*>::iterator bondIter;
265 >    std::vector<Bend*>::iterator bendIter;
266 >    std::vector<Torsion*>::iterator torsionIter;
267 >    Bond* bond;
268 >    Bend* bend;
269 >    Torsion* torsion;
270 >    int a;
271 >    int b;
272 >    int c;
273 >    int d;
274 >    
275 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
276 >        a = bond->getAtomA()->getGlobalIndex();
277 >        b = bond->getAtomB()->getGlobalIndex();        
278 >        exclude_.addPair(a, b);
279 >    }
280  
281 +    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
282 +        a = bend->getAtomA()->getGlobalIndex();
283 +        b = bend->getAtomB()->getGlobalIndex();        
284 +        c = bend->getAtomC()->getGlobalIndex();
285  
286 <  // boxLz
287 <  
288 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
289 <  dsq = dx*dx + dy*dy + dz*dz;
247 <  boxL[2] = sqrt( dsq );
248 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
286 >        exclude_.addPair(a, b);
287 >        exclude_.addPair(a, c);
288 >        exclude_.addPair(b, c);        
289 >    }
290  
291 <  //calculate the max cutoff
292 <  maxCutoff =  calcMaxCutOff();
293 <  
294 <  checkCutOffs();
291 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
292 >        a = torsion->getAtomA()->getGlobalIndex();
293 >        b = torsion->getAtomB()->getGlobalIndex();        
294 >        c = torsion->getAtomC()->getGlobalIndex();        
295 >        d = torsion->getAtomD()->getGlobalIndex();        
296  
297 +        exclude_.addPair(a, b);
298 +        exclude_.addPair(a, c);
299 +        exclude_.addPair(a, d);
300 +        exclude_.addPair(b, c);
301 +        exclude_.addPair(b, d);
302 +        exclude_.addPair(c, d);        
303 +    }
304 +
305 +    
306   }
307  
308 + void SimInfo::removeExcludePairs(Molecule* mol) {
309 +    std::vector<Bond*>::iterator bondIter;
310 +    std::vector<Bend*>::iterator bendIter;
311 +    std::vector<Torsion*>::iterator torsionIter;
312 +    Bond* bond;
313 +    Bend* bend;
314 +    Torsion* torsion;
315 +    int a;
316 +    int b;
317 +    int c;
318 +    int d;
319 +    
320 +    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
321 +        a = bond->getAtomA()->getGlobalIndex();
322 +        b = bond->getAtomB()->getGlobalIndex();        
323 +        exclude_.removePair(a, b);
324 +    }
325  
326 < double SimInfo::calcMaxCutOff(){
326 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
327 >        a = bend->getAtomA()->getGlobalIndex();
328 >        b = bend->getAtomB()->getGlobalIndex();        
329 >        c = bend->getAtomC()->getGlobalIndex();
330  
331 <  double ri[3], rj[3], rk[3];
332 <  double rij[3], rjk[3], rki[3];
333 <  double minDist;
331 >        exclude_.removePair(a, b);
332 >        exclude_.removePair(a, c);
333 >        exclude_.removePair(b, c);        
334 >    }
335  
336 <  ri[0] = Hmat[0][0];
337 <  ri[1] = Hmat[1][0];
338 <  ri[2] = Hmat[2][0];
336 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
337 >        a = torsion->getAtomA()->getGlobalIndex();
338 >        b = torsion->getAtomB()->getGlobalIndex();        
339 >        c = torsion->getAtomC()->getGlobalIndex();        
340 >        d = torsion->getAtomD()->getGlobalIndex();        
341  
342 <  rj[0] = Hmat[0][1];
343 <  rj[1] = Hmat[1][1];
344 <  rj[2] = Hmat[2][1];
342 >        exclude_.removePair(a, b);
343 >        exclude_.removePair(a, c);
344 >        exclude_.removePair(a, d);
345 >        exclude_.removePair(b, c);
346 >        exclude_.removePair(b, d);
347 >        exclude_.removePair(c, d);        
348 >    }
349  
350 <  rk[0] = Hmat[0][2];
273 <  rk[1] = Hmat[1][2];
274 <  rk[2] = Hmat[2][2];
275 <    
276 <  crossProduct3(ri, rj, rij);
277 <  distXY = dotProduct3(rk,rij) / norm3(rij);
350 > }
351  
279  crossProduct3(rj,rk, rjk);
280  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
352  
353 <  crossProduct3(rk,ri, rki);
354 <  distZX = dotProduct3(rj,rki) / norm3(rki);
353 > void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
354 >    int curStampId;
355  
356 <  minDist = min(min(distXY, distYZ), distZX);
357 <  return minDist/2;
358 <  
356 >    //index from 0
357 >    curStampId = molStampIds_.size();
358 >
359 >    moleculeStamps_.push_back(molStamp);
360 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
361   }
362  
363 < void SimInfo::wrapVector( double thePos[3] ){
363 > void SimInfo::update() {
364  
365 <  int i;
366 <  double scaled[3];
365 > #ifdef IS_MPI
366 >    setupFortranParallel();
367 > #endif
368  
369 <  if( !orthoRhombic ){
296 <    // calc the scaled coordinates.
297 <  
369 >    setupFortranSim();
370  
371 <    matVecMul3(HmatInv, thePos, scaled);
372 <    
373 <    for(i=0; i<3; i++)
374 <      scaled[i] -= roundMe(scaled[i]);
303 <    
304 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
305 <    
306 <    matVecMul3(Hmat, scaled, thePos);
371 >    calcNdf();
372 >    calcNdfRaw();
373 >    calcNdfTrans();
374 > }
375  
376 <  }
377 <  else{
378 <    // calc the scaled coordinates.
376 > std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
377 >    typename SimInfo::MoleculeIterator mi;
378 >    Molecule* mol;
379 >    typename Molecule::AtomIterator ai;
380 >    Atom* atom;
381 >    std::set<AtomType*> atomTypes;
382 >
383 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
384 >
385 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
386 >            atomTypes.insert(atom->getAtomType());
387 >        }
388 >        
389 >    }
390 >
391 >    return atomTypes;        
392 > }
393 >
394 > void SimInfo::setupSimType() {
395 >    std::set<AtomType*>::iterator i;
396 >    std::set<AtomType*> atomTypes;
397 >    atomTypes = getUniqueAtomTypes();
398      
399 <    for(i=0; i<3; i++)
400 <      scaled[i] = thePos[i]*HmatInv[i][i];
401 <    
402 <    // wrap the scaled coordinates
403 <    
404 <    for(i=0; i<3; i++)
405 <      scaled[i] -= roundMe(scaled[i]);
406 <    
407 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
408 <    
409 <    for(i=0; i<3; i++)
410 <      thePos[i] = scaled[i]*Hmat[i][i];
411 <  }
412 <    
413 < }
399 >    int useLennardJones = 0;
400 >    int useElectrostatic = 0;
401 >    int useEAM = 0;
402 >    int useCharge = 0;
403 >    int useDirectional = 0;
404 >    int useDipole = 0;
405 >    int useGayBerne = 0;
406 >    int useSticky = 0;
407 >    int useShape = 0;
408 >    int useFLARB = 0; //it is not in AtomType yet
409 >    int useDirectionalAtom = 0;    
410 >    int useElectrostatics = 0;
411 >    //usePBC and useRF are from globals
412 >    bool usePBC = globals_->getPBC();
413 >    bool useRF = globals_->getUseRF();
414  
415 +    //loop over all of the atom types
416 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
417 +        useLennardJones |= i->isLennardJones();
418 +        useElectrostatic |= i->isElectrostatic();
419 +        useEAM |= i->isEAM();
420 +        useCharge |= i->isCharge();
421 +        useDirectional |= i->isDirectional();
422 +        useDipole |= i->isDipole();
423 +        useGayBerne |= i->isGayBerne();
424 +        useSticky |= i->isSticky();
425 +        useShape |= i->isShape();
426 +    }
427  
428 < int SimInfo::getNDF(){
429 <  int ndf_local;
428 >    if (useSticky || useDipole || useGayBerne || useShape) {
429 >        useDirectionalAtom = 1;
430 >    }
431  
432 <  ndf_local = 0;
433 <  
334 <  for(int i = 0; i < integrableObjects.size(); i++){
335 <    ndf_local += 3;
336 <    if (integrableObjects[i]->isDirectional()) {
337 <      if (integrableObjects[i]->isLinear())
338 <        ndf_local += 2;
339 <      else
340 <        ndf_local += 3;
432 >    if (useCharge || useDipole) {
433 >        useElectrostatics = 1;
434      }
342  }
435  
436 <  // n_constraints is local, so subtract them on each processor:
436 > #ifdef IS_MPI    
437 >    int temp;
438  
439 <  ndf_local -= n_constraints;
439 >    temp = usePBC;
440 >    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
441  
442 < #ifdef IS_MPI
443 <  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
350 < #else
351 <  ndf = ndf_local;
352 < #endif
442 >    temp = useDirectionalAtom;
443 >    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
444  
445 <  // nZconstraints is global, as are the 3 COM translations for the
446 <  // entire system:
445 >    temp = useLennardJones;
446 >    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
447  
448 <  ndf = ndf - 3 - nZconstraints;
448 >    temp = useElectrostatics;
449 >    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
450  
451 <  return ndf;
452 < }
451 >    temp = useCharge;
452 >    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
453  
454 < int SimInfo::getNDFraw() {
455 <  int ndfRaw_local;
454 >    temp = useDipole;
455 >    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
456  
457 <  // Raw degrees of freedom that we have to set
458 <  ndfRaw_local = 0;
457 >    temp = useSticky;
458 >    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
459  
460 <  for(int i = 0; i < integrableObjects.size(); i++){
461 <    ndfRaw_local += 3;
370 <    if (integrableObjects[i]->isDirectional()) {
371 <       if (integrableObjects[i]->isLinear())
372 <        ndfRaw_local += 2;
373 <      else
374 <        ndfRaw_local += 3;
375 <    }
376 <  }
377 <    
378 < #ifdef IS_MPI
379 <  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
380 < #else
381 <  ndfRaw = ndfRaw_local;
382 < #endif
460 >    temp = useGayBerne;
461 >    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
462  
463 <  return ndfRaw;
464 < }
463 >    temp = useEAM;
464 >    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
465  
466 < int SimInfo::getNDFtranslational() {
467 <  int ndfTrans_local;
466 >    temp = useShape;
467 >    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
468  
469 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
469 >    temp = useFLARB;
470 >    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
471  
472 <
473 < #ifdef IS_MPI
474 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
395 < #else
396 <  ndfTrans = ndfTrans_local;
472 >    temp = useRF;
473 >    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
474 >    
475   #endif
476  
477 <  ndfTrans = ndfTrans - 3 - nZconstraints;
477 >    fInfo_.SIM_uses_PBC = usePBC;    
478 >    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
479 >    fInfo_.SIM_uses_LennardJones = useLennardJones;
480 >    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
481 >    fInfo_.SIM_uses_Charges = useCharge;
482 >    fInfo_.SIM_uses_Dipoles = useDipole;
483 >    fInfo_.SIM_uses_Sticky = useSticky;
484 >    fInfo_.SIM_uses_GayBerne = useGayBerne;
485 >    fInfo_.SIM_uses_EAM = useEAM;
486 >    fInfo_.SIM_uses_Shapes = useShape;
487 >    fInfo_.SIM_uses_FLARB = useFLARB;
488 >    fInfo_.SIM_uses_RF = useRF;
489  
490 <  return ndfTrans;
490 >    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
491 >        fInfo_.dielect = dielectric;
492 >    } else {
493 >        fInfo_.dielect = 0.0;
494 >    }
495 >
496   }
497  
498 < int SimInfo::getTotIntegrableObjects() {
499 <  int nObjs_local;
500 <  int nObjs;
498 > void SimInfo::setupFortranSim() {
499 >    int isError;
500 >    int nExclude;
501 >    std::vector<int> fortranGlobalGroupMembership;
502 >    
503 >    nExclude = exclude_.getSize();
504 >    isError = 0;
505  
506 <  nObjs_local =  integrableObjects.size();
506 >    //globalGroupMembership_ is filled by SimCreator    
507 >    for (int i = 0; i < nGlobalAtoms_; i++) {
508 >        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
509 >    }
510  
511 +    //calculate mass ratio of cutoff group
512 +    std::vector<double> mfact;
513 +    typename SimInfo::MoleculeIterator mi;
514 +    Molecule* mol;
515 +    typename Molecule::CutoffGroupIterator ci;
516 +    CutoffGroup* cg;
517 +    typename Molecule::AtomIterator ai;
518 +    Atom* atom;
519 +    double totalMass;
520  
521 < #ifdef IS_MPI
522 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
523 < #else
524 <  nObjs = nObjs_local;
525 < #endif
521 >    //to avoid memory reallocation, reserve enough space for mfact
522 >    mfact.reserve(getNCutoffGroups());
523 >    
524 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
525 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
526  
527 +            totalMass = cg->getMass();
528 +            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
529 +                        mfact.push_back(atom->getMass()/totalMass);
530 +            }
531  
532 <  return nObjs;
533 < }
532 >        }      
533 >    }
534  
535 < void SimInfo::refreshSim(){
535 >    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
536 >    std::vector<int> identArray;
537  
538 <  simtype fInfo;
539 <  int isError;
540 <  int n_global;
541 <  int* excl;
538 >    //to avoid memory reallocation, reserve enough space identArray
539 >    identArray.reserve(getNAtoms());
540 >    
541 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
542 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
543 >            identArray.push_back(atom->getIdent());
544 >        }
545 >    }    
546  
547 <  fInfo.dielect = 0.0;
547 >    //fill molMembershipArray
548 >    //molMembershipArray is filled by SimCreator    
549 >    std::vector<int> molMembershipArray(nGlobalAtoms_);
550 >    for (int i = 0; i < nGlobalAtoms_; i++) {
551 >        molMembershipArray.push_back(globalMolMembership_[i] + 1);
552 >    }
553 >    
554 >    //setup fortran simulation
555 >    //gloalExcludes and molMembershipArray should go away (They are never used)
556 >    //why the hell fortran need to know molecule?
557 >    //OOPSE = Object-Obfuscated Parallel Simulation Engine
558 >    
559 >    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
560 >                  &nGlobalExcludes, globalExcludes, molMembershipArray,
561 >                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
562  
563 <  if( useDipoles ){
431 <    if( useReactionField )fInfo.dielect = dielectric;
432 <  }
563 >    if( isError ){
564  
565 <  fInfo.SIM_uses_PBC = usePBC;
566 <  //fInfo.SIM_uses_LJ = 0;
567 <  fInfo.SIM_uses_LJ = useLJ;
568 <  fInfo.SIM_uses_sticky = useSticky;
569 <  //fInfo.SIM_uses_sticky = 0;
570 <  fInfo.SIM_uses_charges = useCharges;
440 <  fInfo.SIM_uses_dipoles = useDipoles;
441 <  //fInfo.SIM_uses_dipoles = 0;
442 <  fInfo.SIM_uses_RF = useReactionField;
443 <  //fInfo.SIM_uses_RF = 0;
444 <  fInfo.SIM_uses_GB = useGB;
445 <  fInfo.SIM_uses_EAM = useEAM;
565 >        sprintf( painCave.errMsg,
566 >                 "There was an error setting the simulation information in fortran.\n" );
567 >        painCave.isFatal = 1;
568 >        painCave.severity = OOPSE_ERROR;
569 >        simError();
570 >    }
571  
447  n_exclude = excludes->getSize();
448  excl = excludes->getFortranArray();
449  
572   #ifdef IS_MPI
573 <  n_global = mpiSim->getNAtomsGlobal();
574 < #else
575 <  n_global = n_atoms;
576 < #endif
577 <  
578 <  isError = 0;
579 <  
458 <  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
459 <  //it may not be a good idea to pass the address of first element in vector
460 <  //since c++ standard does not require vector to be stored continuously in meomory
461 <  //Most of the compilers will organize the memory of vector continuously
462 <  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
463 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
464 <                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
465 <
466 <  if( isError ){
467 <    
468 <    sprintf( painCave.errMsg,
469 <             "There was an error setting the simulation information in fortran.\n" );
470 <    painCave.isFatal = 1;
471 <    painCave.severity = OOPSE_ERROR;
472 <    simError();
473 <  }
474 <  
573 >    sprintf( checkPointMsg,
574 >       "succesfully sent the simulation information to fortran.\n");
575 >    MPIcheckPoint();
576 > #endif // is_mpi
577 > }
578 >
579 >
580   #ifdef IS_MPI
581 <  sprintf( checkPointMsg,
582 <           "succesfully sent the simulation information to fortran.\n");
583 <  MPIcheckPoint();
584 < #endif // is_mpi
585 <  
586 <  this->ndf = this->getNDF();
587 <  this->ndfRaw = this->getNDFraw();
588 <  this->ndfTrans = this->getNDFtranslational();
589 < }
581 > void SimInfo::setupFortranParallel() {
582 >    
583 >    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
584 >    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
585 >    std::vector<int> localToGlobalCutoffGroupIndex;
586 >    typename SimInfo::MoleculeIterator mi;
587 >    typename Molecule::AtomIterator ai;
588 >    typename Molecule::CutoffGroupIterator ci;
589 >    Molecule* mol;
590 >    Atom* atom;
591 >    CutoffGroup* cg;
592 >    mpiSimData parallelData;
593 >    int isError;
594  
595 < void SimInfo::setDefaultRcut( double theRcut ){
487 <  
488 <  haveRcut = 1;
489 <  rCut = theRcut;
490 <  rList = rCut + 1.0;
491 <  
492 <  notifyFortranCutOffs( &rCut, &rSw, &rList );
493 < }
595 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
596  
597 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
597 >        //local index(index in DataStorge) of atom is important
598 >        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
599 >            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
600 >        }
601  
602 <  rSw = theRsw;
603 <  setDefaultRcut( theRcut );
604 < }
602 >        //local index of cutoff group is trivial, it only depends on the order of travesing
603 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
604 >            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
605 >        }        
606 >        
607 >    }
608  
609 +    //fill up mpiSimData struct
610 +    parallelData.nMolGlobal = getNGlobalMolecules();
611 +    parallelData.nMolLocal = getNMolecules();
612 +    parallelData.nAtomsGlobal = getNGlobalAtoms();
613 +    parallelData.nAtomsLocal = getNAtoms();
614 +    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
615 +    parallelData.nGroupsLocal = getNCutoffGroups();
616 +    parallelData.myNode = worldRank;
617 +    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
618  
619 < void SimInfo::checkCutOffs( void ){
620 <  
621 <  if( boxIsInit ){
622 <    
506 <    //we need to check cutOffs against the box
507 <    
508 <    if( rCut > maxCutoff ){
509 <      sprintf( painCave.errMsg,
510 <               "cutoffRadius is too large for the current periodic box.\n"
511 <               "\tCurrent Value of cutoffRadius = %G at time %G\n "
512 <               "\tThis is larger than half of at least one of the\n"
513 <               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
514 <               "\n"
515 <               "\t[ %G %G %G ]\n"
516 <               "\t[ %G %G %G ]\n"
517 <               "\t[ %G %G %G ]\n",
518 <               rCut, currentTime,
519 <               Hmat[0][0], Hmat[0][1], Hmat[0][2],
520 <               Hmat[1][0], Hmat[1][1], Hmat[1][2],
521 <               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
522 <      painCave.severity = OOPSE_ERROR;
523 <      painCave.isFatal = 1;
524 <      simError();
525 <    }    
526 <  } else {
527 <    // initialize this stuff before using it, OK?
528 <    sprintf( painCave.errMsg,
529 <             "Trying to check cutoffs without a box.\n"
530 <             "\tOOPSE should have better programmers than that.\n" );
531 <    painCave.severity = OOPSE_ERROR;
532 <    painCave.isFatal = 1;
533 <    simError();      
534 <  }
535 <  
536 < }
619 >    //pass mpiSimData struct and index arrays to fortran
620 >    setFsimParallel(parallelData, &(parallelData->nAtomsLocal),
621 >                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
622 >                    &localToGlobalCutoffGroupIndex[0], &isError);
623  
624 < void SimInfo::addProperty(GenericData* prop){
624 >    if (isError) {
625 >        sprintf(painCave.errMsg,
626 >                "mpiRefresh errror: fortran didn't like something we gave it.\n");
627 >        painCave.isFatal = 1;
628 >        simError();
629 >    }
630  
631 <  map<string, GenericData*>::iterator result;
632 <  result = properties.find(prop->getID());
542 <  
543 <  //we can't simply use  properties[prop->getID()] = prop,
544 <  //it will cause memory leak if we already contain a propery which has the same name of prop
545 <  
546 <  if(result != properties.end()){
547 <    
548 <    delete (*result).second;
549 <    (*result).second = prop;
550 <      
551 <  }
552 <  else{
631 >    sprintf(checkPointMsg, " mpiRefresh successful.\n");
632 >    MPIcheckPoint();
633  
554    properties[prop->getID()] = prop;
634  
556  }
557    
635   }
636  
637 < GenericData* SimInfo::getProperty(const string& propName){
561 <
562 <  map<string, GenericData*>::iterator result;
563 <  
564 <  //string lowerCaseName = ();
565 <  
566 <  result = properties.find(propName);
567 <  
568 <  if(result != properties.end())
569 <    return (*result).second;  
570 <  else  
571 <    return NULL;  
572 < }
637 > #endif
638  
639 + double SimInfo::calcMaxCutoffRadius() {
640  
575 void SimInfo::getFortranGroupArrays(SimInfo* info,
576                                    vector<int>& FglobalGroupMembership,
577                                    vector<double>& mfact){
578  
579  Molecule* myMols;
580  Atom** myAtoms;
581  int numAtom;
582  double mtot;
583  int numMol;
584  int numCutoffGroups;
585  CutoffGroup* myCutoffGroup;
586  vector<CutoffGroup*>::iterator iterCutoff;
587  Atom* cutoffAtom;
588  vector<Atom*>::iterator iterAtom;
589  int atomIndex;
590  double totalMass;
591  
592  mfact.clear();
593  FglobalGroupMembership.clear();
594  
641  
642 <  // Fix the silly fortran indexing problem
642 >    std::vector<AtomType*> atomTypes;
643 >    std::vector<AtomType*>::iterator i;
644 >    std::vector<double> cutoffRadius;
645 >
646 >    //get the unique atom types
647 >    atomTypes = getUniqueAtomTypes();
648 >
649 >    //query the max cutoff radius among these atom types
650 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
651 >        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
652 >    }
653 >
654 >    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
655   #ifdef IS_MPI
656 <  numAtom = mpiSim->getNAtomsGlobal();
599 < #else
600 <  numAtom = n_atoms;
656 >    //pick the max cutoff radius among the processors
657   #endif
602  for (int i = 0; i < numAtom; i++)
603    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
604  
658  
659 <  myMols = info->molecules;
660 <  numMol = info->n_mol;
608 <  for(int i  = 0; i < numMol; i++){
609 <    numCutoffGroups = myMols[i].getNCutoffGroups();
610 <    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
611 <        myCutoffGroup != NULL;
612 <        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
659 >    return maxCutoffRadius;
660 > }
661  
662 <      totalMass = myCutoffGroup->getMass();
662 > void SimInfo::addProperty(GenericData* genData) {
663 >    properties_.addProperty(genData);  
664 > }
665 >
666 > void SimInfo::removeProperty(const std::string& propName) {
667 >    properties_.removeProperty(propName);  
668 > }
669 >
670 > void SimInfo::clearProperties() {
671 >    properties_.clearProperties();
672 > }
673 >
674 > std::vector<std::string> SimInfo::getPropertyNames() {
675 >    return properties_.getPropertyNames();  
676 > }
677        
678 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
679 <          cutoffAtom != NULL;
680 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
619 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
620 <      }  
621 <    }
622 <  }
678 > std::vector<GenericData*> SimInfo::getProperties() {
679 >    return properties_.getProperties();
680 > }
681  
682 + GenericData* SimInfo::getPropertyByName(const std::string& propName) {
683 +    return properties_.getPropertyByName(propName);
684   }
685 +
686 +
687 + std::ostream& operator <<(ostream& o, SimInfo& info) {
688 +
689 +    return o;
690 + }
691 +
692 + }//end namespace oopse

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines