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 1886 by tim, Wed Dec 15 16:13:35 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"
34 < #define __C
35 < #include "fSimulation.h"
36 < #include "simError.h"
33 > #include <algorithm>
34 > #include <set>
35 > #include "brains/SimInfo.hpp"
36 > #include "primitives/Molecule.hpp"
37 > #include "UseTheForce/doForces_interface.h"
38 > #include "UseTheForce/notifyCutoffs_interface.h"
39 > #include "utils/MemoryUtils.hpp"
40 > #include "utils/simError.h"
41  
13 – #include "fortranWrappers.hpp"
14 –
15 – #include "MatVec3.h"
16 –
42   #ifdef IS_MPI
43 < #include "mpiSimulation.hpp"
44 < #endif
43 > #include "UseTheForce/mpiComponentPlan.h"
44 > #include "UseTheForce/DarkSide/simParallel_interface.h"
45 > #endif
46  
47 < inline double roundMe( double x ){
22 <  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
23 < }
24 <          
25 < inline double min( double a, double b ){
26 <  return (a < b ) ? a : b;
27 < }
47 > namespace oopse {
48  
49 < SimInfo* currentInfo;
49 > SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
50 >                                ForceField* ff, Globals* simParams) :
51 >                                forceField_(ff), simParams_(simParams),
52 >                                ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
53 >                                nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
54 >                                nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
55 >                                nAtoms_(0), nBonds_(0),  nBends_(0), nTorsions_(0), nRigidBodies_(0),
56 >                                nIntegrableObjects_(0),  nCutoffGroups_(0), nConstraints_(0),
57 >                                sman_(NULL), fortranInitialized_(false) {
58  
59 < SimInfo::SimInfo(){
59 >            
60 >    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
61 >    MoleculeStamp* molStamp;
62 >    int nMolWithSameStamp;
63 >    int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
64 >    int nGroups = 0;          //total cutoff groups defined in meta-data file
65 >    CutoffGroupStamp* cgStamp;    
66 >    RigidBodyStamp* rbStamp;
67 >    int nRigidAtoms = 0;
68 >    
69 >    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
70 >        molStamp = i->first;
71 >        nMolWithSameStamp = i->second;
72 >        
73 >        addMoleculeStamp(molStamp, nMolWithSameStamp);
74  
75 <  n_constraints = 0;
76 <  nZconstraints = 0;
35 <  n_oriented = 0;
36 <  n_dipoles = 0;
37 <  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;
75 >        //calculate atoms in molecules
76 >        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
77  
47 –  haveRcut = 0;
48 –  haveRsw = 0;
49 –  boxIsInit = 0;
50 –  
51 –  resetTime = 1e99;
78  
79 <  orthoRhombic = 0;
80 <  orthoTolerance = 1E-6;
81 <  useInitXSstate = true;
79 >        //calculate atoms in cutoff groups
80 >        int nAtomsInGroups = 0;
81 >        int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
82 >        
83 >        for (int j=0; j < nCutoffGroupsInStamp; j++) {
84 >            cgStamp = molStamp->getCutoffGroup(j);
85 >            nAtomsInGroups += cgStamp->getNMembers();
86 >        }
87  
88 <  usePBC = 0;
89 <  useLJ = 0;
59 <  useSticky = 0;
60 <  useCharges = 0;
61 <  useDipoles = 0;
62 <  useReactionField = 0;
63 <  useGB = 0;
64 <  useEAM = 0;
65 <  useSolidThermInt = 0;
66 <  useLiquidThermInt = 0;
88 >        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
89 >        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
90  
91 <  haveCutoffGroups = false;
91 >        //calculate atoms in rigid bodies
92 >        int nAtomsInRigidBodies = 0;
93 >        int nRigidBodiesInStamp = molStamp->getNCutoffGroups();
94 >        
95 >        for (int j=0; j < nRigidBodiesInStamp; j++) {
96 >            rbStamp = molStamp->getRigidBody(j);
97 >            nAtomsInRigidBodies += rbStamp->getNMembers();
98 >        }
99  
100 <  excludes = Exclude::Instance();
100 >        nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
101 >        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
102 >        
103 >    }
104  
105 <  myConfiguration = new SimState();
105 >    //every free atom (atom does not belong to cutoff groups) is a cutoff group
106 >    //therefore the total number of cutoff groups in the system is equal to
107 >    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
108 >    //file plus the number of cutoff groups defined in meta-data file
109 >    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
110  
111 <  has_minimizer = false;
112 <  the_minimizer =NULL;
111 >    //every free atom (atom does not belong to rigid bodies) is an integrable object
112 >    //therefore the total number of  integrable objects in the system is equal to
113 >    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
114 >    //file plus the number of  rigid bodies defined in meta-data file
115 >    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
116  
117 <  ngroup = 0;
117 >    nGlobalMols_ = molStampIds_.size();
118  
119 <  wrapMeSimInfo( this );
119 > #ifdef IS_MPI    
120 >    molToProcMap_.resize(nGlobalMols_);
121 > #endif
122 >    
123   }
124  
125 + SimInfo::~SimInfo() {
126 +    //MemoryUtils::deleteVectorOfPointer(molecules_);
127  
128 < SimInfo::~SimInfo(){
128 >    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
129 >    
130 >    delete sman_;
131 >    delete simParams_;
132 >    delete forceField_;
133  
85 –  delete myConfiguration;
86 –
87 –  map<string, GenericData*>::iterator i;
88 –  
89 –  for(i = properties.begin(); i != properties.end(); i++)
90 –    delete (*i).second;
91 –
134   }
135  
94 – void SimInfo::setBox(double newBox[3]) {
95 –  
96 –  int i, j;
97 –  double tempMat[3][3];
136  
137 <  for(i=0; i<3; i++)
138 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
137 > bool SimInfo::addMolecule(Molecule* mol) {
138 >    MoleculeIterator i;
139  
140 <  tempMat[0][0] = newBox[0];
141 <  tempMat[1][1] = newBox[1];
104 <  tempMat[2][2] = newBox[2];
140 >    i = molecules_.find(mol->getGlobalIndex());
141 >    if (i == molecules_.end() ) {
142  
143 <  setBoxM( tempMat );
143 >        molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
144 >        
145 >        nAtoms_ += mol->getNAtoms();
146 >        nBonds_ += mol->getNBonds();
147 >        nBends_ += mol->getNBends();
148 >        nTorsions_ += mol->getNTorsions();
149 >        nRigidBodies_ += mol->getNRigidBodies();
150 >        nIntegrableObjects_ += mol->getNIntegrableObjects();
151 >        nCutoffGroups_ += mol->getNCutoffGroups();
152 >        nConstraints_ += mol->getNConstraints();
153  
154 +        addExcludePairs(mol);
155 +        
156 +        return true;
157 +    } else {
158 +        return false;
159 +    }
160   }
161  
162 < void SimInfo::setBoxM( double theBox[3][3] ){
163 <  
164 <  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);
162 > bool SimInfo::removeMolecule(Molecule* mol) {
163 >    MoleculeIterator i;
164 >    i = molecules_.find(mol->getGlobalIndex());
165  
166 <  if( !boxIsInit ) boxIsInit = 1;
166 >    if (i != molecules_.end() ) {
167  
168 <  for(i=0; i < 3; i++)
169 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
170 <  
171 <  calcBoxL();
172 <  calcHmatInv();
168 >        assert(mol == i->second);
169 >        
170 >        nAtoms_ -= mol->getNAtoms();
171 >        nBonds_ -= mol->getNBonds();
172 >        nBends_ -= mol->getNBends();
173 >        nTorsions_ -= mol->getNTorsions();
174 >        nRigidBodies_ -= mol->getNRigidBodies();
175 >        nIntegrableObjects_ -= mol->getNIntegrableObjects();
176 >        nCutoffGroups_ -= mol->getNCutoffGroups();
177 >        nConstraints_ -= mol->getNConstraints();
178  
179 <  for(i=0; i < 3; i++) {
180 <    for (j=0; j < 3; j++) {
181 <      FortranHmat[3*j + i] = Hmat[i][j];
182 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
179 >        removeExcludePairs(mol);
180 >        molecules_.erase(mol->getGlobalIndex());
181 >
182 >        delete mol;
183 >        
184 >        return true;
185 >    } else {
186 >        return false;
187      }
133 –  }
188  
135 –  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
136 –
137 – }
138 –
189  
190 < void SimInfo::getBoxM (double theBox[3][3]) {
190 > }    
191  
192 <  int i, j;
193 <  for(i=0; i<3; i++)
194 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
192 >        
193 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
194 >    i = molecules_.begin();
195 >    return i == molecules_.end() ? NULL : i->second;
196 > }    
197 >
198 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
199 >    ++i;
200 >    return i == molecules_.end() ? NULL : i->second;    
201   }
202  
203  
204 < void SimInfo::scaleBox(double scale) {
205 <  double theBox[3][3];
206 <  int i, j;
204 > void SimInfo::calcNdf() {
205 >    int ndf_local;
206 >    MoleculeIterator i;
207 >    std::vector<StuntDouble*>::iterator j;
208 >    Molecule* mol;
209 >    StuntDouble* integrableObject;
210  
211 <  // cerr << "Scaling box by " << scale << "\n";
211 >    ndf_local = 0;
212 >    
213 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
214 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
215 >               integrableObject = mol->nextIntegrableObject(j)) {
216  
217 <  for(i=0; i<3; i++)
155 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
217 >            ndf_local += 3;
218  
219 <  setBoxM(theBox);
219 >            if (integrableObject->isDirectional()) {
220 >                if (integrableObject->isLinear()) {
221 >                    ndf_local += 2;
222 >                } else {
223 >                    ndf_local += 3;
224 >                }
225 >            }
226 >            
227 >        }//end for (integrableObject)
228 >    }// end for (mol)
229 >    
230 >    // n_constraints is local, so subtract them on each processor
231 >    ndf_local -= nConstraints_;
232  
233 + #ifdef IS_MPI
234 +    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
235 + #else
236 +    ndf_ = ndf_local;
237 + #endif
238 +
239 +    // nZconstraints_ is global, as are the 3 COM translations for the
240 +    // entire system:
241 +    ndf_ = ndf_ - 3 - nZconstraint_;
242 +
243   }
244  
245 < void SimInfo::calcHmatInv( void ) {
246 <  
163 <  int oldOrtho;
164 <  int i,j;
165 <  double smallDiag;
166 <  double tol;
167 <  double sanity[3][3];
245 > void SimInfo::calcNdfRaw() {
246 >    int ndfRaw_local;
247  
248 <  invertMat3( Hmat, HmatInv );
248 >    MoleculeIterator i;
249 >    std::vector<StuntDouble*>::iterator j;
250 >    Molecule* mol;
251 >    StuntDouble* integrableObject;
252  
253 <  // check to see if Hmat is orthorhombic
254 <  
255 <  oldOrtho = orthoRhombic;
253 >    // Raw degrees of freedom that we have to set
254 >    ndfRaw_local = 0;
255 >    
256 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
257 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
258 >               integrableObject = mol->nextIntegrableObject(j)) {
259  
260 <  smallDiag = fabs(Hmat[0][0]);
176 <  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;
260 >            ndfRaw_local += 3;
261  
262 <  orthoRhombic = 1;
263 <  
264 <  for (i = 0; i < 3; i++ ) {
265 <    for (j = 0 ; j < 3; j++) {
266 <      if (i != j) {
267 <        if (orthoRhombic) {
268 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
269 <        }        
270 <      }
262 >            if (integrableObject->isDirectional()) {
263 >                if (integrableObject->isLinear()) {
264 >                    ndfRaw_local += 2;
265 >                } else {
266 >                    ndfRaw_local += 3;
267 >                }
268 >            }
269 >            
270 >        }
271      }
190 –  }
191 –
192 –  if( oldOrtho != orthoRhombic ){
272      
273 <    if( orthoRhombic ) {
274 <      sprintf( painCave.errMsg,
275 <               "OOPSE is switching from the default Non-Orthorhombic\n"
276 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
277 <               "\tThis is usually a good thing, but if you wan't the\n"
199 <               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
200 <               "\tvariable ( currently set to %G ) smaller.\n",
201 <               orthoTolerance);
202 <      painCave.severity = OOPSE_INFO;
203 <      simError();
204 <    }
205 <    else {
206 <      sprintf( painCave.errMsg,
207 <               "OOPSE is switching from the faster Orthorhombic to the more\n"
208 <               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
209 <               "\tThis is usually because the box has deformed under\n"
210 <               "\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 <  }
273 > #ifdef IS_MPI
274 >    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
275 > #else
276 >    ndfRaw_ = ndfRaw_local;
277 > #endif
278   }
279  
280 < void SimInfo::calcBoxL( void ){
280 > void SimInfo::calcNdfTrans() {
281 >    int ndfTrans_local;
282  
283 <  double dx, dy, dz, dsq;
283 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
284  
224 –  // boxVol = Determinant of Hmat
285  
286 <  boxVol = matDet3( Hmat );
286 > #ifdef IS_MPI
287 >    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
288 > #else
289 >    ndfTrans_ = ndfTrans_local;
290 > #endif
291  
292 <  // boxLx
293 <  
294 <  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];
292 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
293 >
294 > }
295  
296 <  // boxLy
297 <  
298 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
299 <  dsq = dx*dx + dy*dy + dz*dz;
300 <  boxL[1] = sqrt( dsq );
301 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
296 > void SimInfo::addExcludePairs(Molecule* mol) {
297 >    std::vector<Bond*>::iterator bondIter;
298 >    std::vector<Bend*>::iterator bendIter;
299 >    std::vector<Torsion*>::iterator torsionIter;
300 >    Bond* bond;
301 >    Bend* bend;
302 >    Torsion* torsion;
303 >    int a;
304 >    int b;
305 >    int c;
306 >    int d;
307 >    
308 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
309 >        a = bond->getAtomA()->getGlobalIndex();
310 >        b = bond->getAtomB()->getGlobalIndex();        
311 >        exclude_.addPair(a, b);
312 >    }
313  
314 +    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
315 +        a = bend->getAtomA()->getGlobalIndex();
316 +        b = bend->getAtomB()->getGlobalIndex();        
317 +        c = bend->getAtomC()->getGlobalIndex();
318  
319 <  // boxLz
320 <  
321 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
322 <  dsq = dx*dx + dy*dy + dz*dz;
247 <  boxL[2] = sqrt( dsq );
248 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
319 >        exclude_.addPair(a, b);
320 >        exclude_.addPair(a, c);
321 >        exclude_.addPair(b, c);        
322 >    }
323  
324 <  //calculate the max cutoff
325 <  maxCutoff =  calcMaxCutOff();
326 <  
327 <  checkCutOffs();
324 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
325 >        a = torsion->getAtomA()->getGlobalIndex();
326 >        b = torsion->getAtomB()->getGlobalIndex();        
327 >        c = torsion->getAtomC()->getGlobalIndex();        
328 >        d = torsion->getAtomD()->getGlobalIndex();        
329  
330 +        exclude_.addPair(a, b);
331 +        exclude_.addPair(a, c);
332 +        exclude_.addPair(a, d);
333 +        exclude_.addPair(b, c);
334 +        exclude_.addPair(b, d);
335 +        exclude_.addPair(c, d);        
336 +    }
337 +
338 +    
339   }
340  
341 + void SimInfo::removeExcludePairs(Molecule* mol) {
342 +    std::vector<Bond*>::iterator bondIter;
343 +    std::vector<Bend*>::iterator bendIter;
344 +    std::vector<Torsion*>::iterator torsionIter;
345 +    Bond* bond;
346 +    Bend* bend;
347 +    Torsion* torsion;
348 +    int a;
349 +    int b;
350 +    int c;
351 +    int d;
352 +    
353 +    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
354 +        a = bond->getAtomA()->getGlobalIndex();
355 +        b = bond->getAtomB()->getGlobalIndex();        
356 +        exclude_.removePair(a, b);
357 +    }
358  
359 < double SimInfo::calcMaxCutOff(){
359 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
360 >        a = bend->getAtomA()->getGlobalIndex();
361 >        b = bend->getAtomB()->getGlobalIndex();        
362 >        c = bend->getAtomC()->getGlobalIndex();
363  
364 <  double ri[3], rj[3], rk[3];
365 <  double rij[3], rjk[3], rki[3];
366 <  double minDist;
364 >        exclude_.removePair(a, b);
365 >        exclude_.removePair(a, c);
366 >        exclude_.removePair(b, c);        
367 >    }
368  
369 <  ri[0] = Hmat[0][0];
370 <  ri[1] = Hmat[1][0];
371 <  ri[2] = Hmat[2][0];
369 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
370 >        a = torsion->getAtomA()->getGlobalIndex();
371 >        b = torsion->getAtomB()->getGlobalIndex();        
372 >        c = torsion->getAtomC()->getGlobalIndex();        
373 >        d = torsion->getAtomD()->getGlobalIndex();        
374  
375 <  rj[0] = Hmat[0][1];
376 <  rj[1] = Hmat[1][1];
377 <  rj[2] = Hmat[2][1];
375 >        exclude_.removePair(a, b);
376 >        exclude_.removePair(a, c);
377 >        exclude_.removePair(a, d);
378 >        exclude_.removePair(b, c);
379 >        exclude_.removePair(b, d);
380 >        exclude_.removePair(c, d);        
381 >    }
382  
383 <  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);
383 > }
384  
279 –  crossProduct3(rj,rk, rjk);
280 –  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
385  
386 <  crossProduct3(rk,ri, rki);
387 <  distZX = dotProduct3(rj,rki) / norm3(rki);
386 > void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
387 >    int curStampId;
388  
389 <  minDist = min(min(distXY, distYZ), distZX);
390 <  return minDist/2;
391 <  
389 >    //index from 0
390 >    curStampId = molStampIds_.size();
391 >
392 >    moleculeStamps_.push_back(molStamp);
393 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
394   }
395  
396 < void SimInfo::wrapVector( double thePos[3] ){
396 > void SimInfo::update() {
397  
398 <  int i;
293 <  double scaled[3];
398 >    setupSimType();
399  
400 <  if( !orthoRhombic ){
401 <    // calc the scaled coordinates.
402 <  
400 > #ifdef IS_MPI
401 >    setupFortranParallel();
402 > #endif
403  
404 <    matVecMul3(HmatInv, thePos, scaled);
300 <    
301 <    for(i=0; i<3; i++)
302 <      scaled[i] -= roundMe(scaled[i]);
303 <    
304 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
305 <    
306 <    matVecMul3(Hmat, scaled, thePos);
404 >    setupFortranSim();
405  
406 <  }
407 <  else{
408 <    // calc the scaled coordinates.
406 >    //setup fortran force field
407 >    /** @deprecate */    
408 >    int isError = 0;
409 >    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
410 >    if(isError){
411 >        sprintf( painCave.errMsg,
412 >         "ForceField error: There was an error initializing the forceField in fortran.\n" );
413 >        painCave.isFatal = 1;
414 >        simError();
415 >    }
416 >  
417      
418 <    for(i=0; i<3; i++)
419 <      scaled[i] = thePos[i]*HmatInv[i][i];
420 <    
421 <    // wrap the scaled coordinates
422 <    
423 <    for(i=0; i<3; i++)
424 <      scaled[i] -= roundMe(scaled[i]);
319 <    
320 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
321 <    
322 <    for(i=0; i<3; i++)
323 <      thePos[i] = scaled[i]*Hmat[i][i];
324 <  }
325 <    
418 >    setupCutoff();
419 >
420 >    calcNdf();
421 >    calcNdfRaw();
422 >    calcNdfTrans();
423 >
424 >    fortranInitialized_ = true;
425   }
426  
427 + std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
428 +    SimInfo::MoleculeIterator mi;
429 +    Molecule* mol;
430 +    Molecule::AtomIterator ai;
431 +    Atom* atom;
432 +    std::set<AtomType*> atomTypes;
433  
434 < int SimInfo::getNDF(){
330 <  int ndf_local;
434 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
435  
436 <  ndf_local = 0;
437 <  
438 <  for(int i = 0; i < integrableObjects.size(); i++){
439 <    ndf_local += 3;
336 <    if (integrableObjects[i]->isDirectional()) {
337 <      if (integrableObjects[i]->isLinear())
338 <        ndf_local += 2;
339 <      else
340 <        ndf_local += 3;
436 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
437 >            atomTypes.insert(atom->getAtomType());
438 >        }
439 >        
440      }
342 –  }
441  
442 <  // n_constraints is local, so subtract them on each processor:
442 >    return atomTypes;        
443 > }
444  
445 <  ndf_local -= n_constraints;
445 > void SimInfo::setupSimType() {
446 >    std::set<AtomType*>::iterator i;
447 >    std::set<AtomType*> atomTypes;
448 >    atomTypes = getUniqueAtomTypes();
449 >    
450 >    int useLennardJones = 0;
451 >    int useElectrostatic = 0;
452 >    int useEAM = 0;
453 >    int useCharge = 0;
454 >    int useDirectional = 0;
455 >    int useDipole = 0;
456 >    int useGayBerne = 0;
457 >    int useSticky = 0;
458 >    int useShape = 0;
459 >    int useFLARB = 0; //it is not in AtomType yet
460 >    int useDirectionalAtom = 0;    
461 >    int useElectrostatics = 0;
462 >    //usePBC and useRF are from simParams
463 >    int usePBC = simParams_->getPBC();
464 >    int useRF = simParams_->getUseRF();
465  
466 < #ifdef IS_MPI
467 <  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
468 < #else
469 <  ndf = ndf_local;
470 < #endif
466 >    //loop over all of the atom types
467 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
468 >        useLennardJones |= (*i)->isLennardJones();
469 >        useElectrostatic |= (*i)->isElectrostatic();
470 >        useEAM |= (*i)->isEAM();
471 >        useCharge |= (*i)->isCharge();
472 >        useDirectional |= (*i)->isDirectional();
473 >        useDipole |= (*i)->isDipole();
474 >        useGayBerne |= (*i)->isGayBerne();
475 >        useSticky |= (*i)->isSticky();
476 >        useShape |= (*i)->isShape();
477 >    }
478  
479 <  // nZconstraints is global, as are the 3 COM translations for the
480 <  // entire system:
479 >    if (useSticky || useDipole || useGayBerne || useShape) {
480 >        useDirectionalAtom = 1;
481 >    }
482  
483 <  ndf = ndf - 3 - nZconstraints;
483 >    if (useCharge || useDipole) {
484 >        useElectrostatics = 1;
485 >    }
486  
487 <  return ndf;
488 < }
487 > #ifdef IS_MPI    
488 >    int temp;
489  
490 < int SimInfo::getNDFraw() {
491 <  int ndfRaw_local;
490 >    temp = usePBC;
491 >    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
492  
493 <  // Raw degrees of freedom that we have to set
494 <  ndfRaw_local = 0;
493 >    temp = useDirectionalAtom;
494 >    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
495  
496 <  for(int i = 0; i < integrableObjects.size(); i++){
497 <    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
496 >    temp = useLennardJones;
497 >    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
498  
499 <  return ndfRaw;
500 < }
499 >    temp = useElectrostatics;
500 >    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
501  
502 < int SimInfo::getNDFtranslational() {
503 <  int ndfTrans_local;
502 >    temp = useCharge;
503 >    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
504  
505 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
505 >    temp = useDipole;
506 >    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
507  
508 +    temp = useSticky;
509 +    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
510  
511 < #ifdef IS_MPI
512 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
395 < #else
396 <  ndfTrans = ndfTrans_local;
397 < #endif
511 >    temp = useGayBerne;
512 >    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
513  
514 <  ndfTrans = ndfTrans - 3 - nZconstraints;
514 >    temp = useEAM;
515 >    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
516  
517 <  return ndfTrans;
518 < }
517 >    temp = useShape;
518 >    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
519  
520 < int SimInfo::getTotIntegrableObjects() {
521 <  int nObjs_local;
406 <  int nObjs;
520 >    temp = useFLARB;
521 >    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
522  
523 <  nObjs_local =  integrableObjects.size();
523 >    temp = useRF;
524 >    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
525 >    
526 > #endif
527  
528 +    fInfo_.SIM_uses_PBC = usePBC;    
529 +    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
530 +    fInfo_.SIM_uses_LennardJones = useLennardJones;
531 +    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
532 +    fInfo_.SIM_uses_Charges = useCharge;
533 +    fInfo_.SIM_uses_Dipoles = useDipole;
534 +    fInfo_.SIM_uses_Sticky = useSticky;
535 +    fInfo_.SIM_uses_GayBerne = useGayBerne;
536 +    fInfo_.SIM_uses_EAM = useEAM;
537 +    fInfo_.SIM_uses_Shapes = useShape;
538 +    fInfo_.SIM_uses_FLARB = useFLARB;
539 +    fInfo_.SIM_uses_RF = useRF;
540  
541 < #ifdef IS_MPI
412 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
413 < #else
414 <  nObjs = nObjs_local;
415 < #endif
541 >    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
542  
543 +        if (simParams_->haveDielectric()) {
544 +            fInfo_.dielect = simParams_->getDielectric();
545 +        } else {
546 +            sprintf(painCave.errMsg,
547 +                    "SimSetup Error: No Dielectric constant was set.\n"
548 +                    "\tYou are trying to use Reaction Field without"
549 +                    "\tsetting a dielectric constant!\n");
550 +            painCave.isFatal = 1;
551 +            simError();
552 +        }
553 +        
554 +    } else {
555 +        fInfo_.dielect = 0.0;
556 +    }
557  
418 –  return nObjs;
558   }
559  
560 < void SimInfo::refreshSim(){
560 > void SimInfo::setupFortranSim() {
561 >    int isError;
562 >    int nExclude;
563 >    std::vector<int> fortranGlobalGroupMembership;
564 >    
565 >    nExclude = exclude_.getSize();
566 >    isError = 0;
567  
568 <  simtype fInfo;
569 <  int isError;
570 <  int n_global;
571 <  int* excl;
568 >    //globalGroupMembership_ is filled by SimCreator    
569 >    for (int i = 0; i < nGlobalAtoms_; i++) {
570 >        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
571 >    }
572  
573 <  fInfo.dielect = 0.0;
573 >    //calculate mass ratio of cutoff group
574 >    std::vector<double> mfact;
575 >    SimInfo::MoleculeIterator mi;
576 >    Molecule* mol;
577 >    Molecule::CutoffGroupIterator ci;
578 >    CutoffGroup* cg;
579 >    Molecule::AtomIterator ai;
580 >    Atom* atom;
581 >    double totalMass;
582  
583 <  if( useDipoles ){
584 <    if( useReactionField )fInfo.dielect = dielectric;
585 <  }
583 >    //to avoid memory reallocation, reserve enough space for mfact
584 >    mfact.reserve(getNCutoffGroups());
585 >    
586 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
587 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
588  
589 <  fInfo.SIM_uses_PBC = usePBC;
590 <  //fInfo.SIM_uses_LJ = 0;
591 <  fInfo.SIM_uses_LJ = useLJ;
592 <  fInfo.SIM_uses_sticky = useSticky;
438 <  //fInfo.SIM_uses_sticky = 0;
439 <  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;
589 >            totalMass = cg->getMass();
590 >            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
591 >                        mfact.push_back(atom->getMass()/totalMass);
592 >            }
593  
594 <  n_exclude = excludes->getSize();
595 <  excl = excludes->getFortranArray();
449 <  
450 < #ifdef IS_MPI
451 <  n_global = mpiSim->getNAtomsGlobal();
452 < #else
453 <  n_global = n_atoms;
454 < #endif
455 <  
456 <  isError = 0;
457 <  
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);
594 >        }      
595 >    }
596  
597 <  if( isError ){
597 >    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
598 >    std::vector<int> identArray;
599 >
600 >    //to avoid memory reallocation, reserve enough space identArray
601 >    identArray.reserve(getNAtoms());
602      
603 <    sprintf( painCave.errMsg,
604 <             "There was an error setting the simulation information in fortran.\n" );
605 <    painCave.isFatal = 1;
606 <    painCave.severity = OOPSE_ERROR;
607 <    simError();
608 <  }
609 <  
603 >    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
604 >        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
605 >            identArray.push_back(atom->getIdent());
606 >        }
607 >    }    
608 >
609 >    //fill molMembershipArray
610 >    //molMembershipArray is filled by SimCreator    
611 >    std::vector<int> molMembershipArray(nGlobalAtoms_);
612 >    for (int i = 0; i < nGlobalAtoms_; i++) {
613 >        molMembershipArray[i] = globalMolMembership_[i] + 1;
614 >    }
615 >    
616 >    //setup fortran simulation
617 >    //gloalExcludes and molMembershipArray should go away (They are never used)
618 >    //why the hell fortran need to know molecule?
619 >    //OOPSE = Object-Obfuscated Parallel Simulation Engine
620 >    int nGlobalExcludes = 0;
621 >    int* globalExcludes = NULL;
622 >    int* excludeList = exclude_.getExcludeList();
623 >    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
624 >                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
625 >                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
626 >
627 >    if( isError ){
628 >
629 >        sprintf( painCave.errMsg,
630 >                 "There was an error setting the simulation information in fortran.\n" );
631 >        painCave.isFatal = 1;
632 >        painCave.severity = OOPSE_ERROR;
633 >        simError();
634 >    }
635 >
636   #ifdef IS_MPI
637 <  sprintf( checkPointMsg,
638 <           "succesfully sent the simulation information to fortran.\n");
639 <  MPIcheckPoint();
637 >    sprintf( checkPointMsg,
638 >       "succesfully sent the simulation information to fortran.\n");
639 >    MPIcheckPoint();
640   #endif // is_mpi
480 –  
481 –  this->ndf = this->getNDF();
482 –  this->ndfRaw = this->getNDFraw();
483 –  this->ndfTrans = this->getNDFtranslational();
641   }
642  
486 – void SimInfo::setDefaultRcut( double theRcut ){
487 –  
488 –  haveRcut = 1;
489 –  rCut = theRcut;
490 –  rList = rCut + 1.0;
491 –  
492 –  notifyFortranCutOffs( &rCut, &rSw, &rList );
493 – }
643  
644 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
644 > #ifdef IS_MPI
645 > void SimInfo::setupFortranParallel() {
646 >    
647 >    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
648 >    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
649 >    std::vector<int> localToGlobalCutoffGroupIndex;
650 >    SimInfo::MoleculeIterator mi;
651 >    Molecule::AtomIterator ai;
652 >    Molecule::CutoffGroupIterator ci;
653 >    Molecule* mol;
654 >    Atom* atom;
655 >    CutoffGroup* cg;
656 >    mpiSimData parallelData;
657 >    int isError;
658  
659 <  rSw = theRsw;
498 <  setDefaultRcut( theRcut );
499 < }
659 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
660  
661 +        //local index(index in DataStorge) of atom is important
662 +        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
663 +            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
664 +        }
665  
666 < void SimInfo::checkCutOffs( void ){
667 <  
668 <  if( boxIsInit ){
669 <    
670 <    //we need to check cutOffs against the box
671 <    
672 <    if( rCut > maxCutoff ){
673 <      sprintf( painCave.errMsg,
674 <               "cutoffRadius is too large for the current periodic box.\n"
675 <               "\tCurrent Value of cutoffRadius = %G at time %G\n "
676 <               "\tThis is larger than half of at least one of the\n"
677 <               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
678 <               "\n"
679 <               "\t[ %G %G %G ]\n"
680 <               "\t[ %G %G %G ]\n"
681 <               "\t[ %G %G %G ]\n",
682 <               rCut, currentTime,
683 <               Hmat[0][0], Hmat[0][1], Hmat[0][2],
684 <               Hmat[1][0], Hmat[1][1], Hmat[1][2],
685 <               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
686 <      painCave.severity = OOPSE_ERROR;
687 <      painCave.isFatal = 1;
688 <      simError();
689 <    }    
690 <  } else {
691 <    // initialize this stuff before using it, OK?
692 <    sprintf( painCave.errMsg,
693 <             "Trying to check cutoffs without a box.\n"
694 <             "\tOOPSE should have better programmers than that.\n" );
695 <    painCave.severity = OOPSE_ERROR;
696 <    painCave.isFatal = 1;
697 <    simError();      
698 <  }
535 <  
666 >        //local index of cutoff group is trivial, it only depends on the order of travesing
667 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
668 >            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
669 >        }        
670 >        
671 >    }
672 >
673 >    //fill up mpiSimData struct
674 >    parallelData.nMolGlobal = getNGlobalMolecules();
675 >    parallelData.nMolLocal = getNMolecules();
676 >    parallelData.nAtomsGlobal = getNGlobalAtoms();
677 >    parallelData.nAtomsLocal = getNAtoms();
678 >    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
679 >    parallelData.nGroupsLocal = getNCutoffGroups();
680 >    parallelData.myNode = worldRank;
681 >    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
682 >
683 >    //pass mpiSimData struct and index arrays to fortran
684 >    setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
685 >                    &localToGlobalAtomIndex[0],  &(parallelData.nGroupsLocal),
686 >                    &localToGlobalCutoffGroupIndex[0], &isError);
687 >
688 >    if (isError) {
689 >        sprintf(painCave.errMsg,
690 >                "mpiRefresh errror: fortran didn't like something we gave it.\n");
691 >        painCave.isFatal = 1;
692 >        simError();
693 >    }
694 >
695 >    sprintf(checkPointMsg, " mpiRefresh successful.\n");
696 >    MPIcheckPoint();
697 >
698 >
699   }
700  
701 < void SimInfo::addProperty(GenericData* prop){
701 > #endif
702  
703 <  map<string, GenericData*>::iterator result;
704 <  result = properties.find(prop->getID());
705 <  
706 <  //we can't simply use  properties[prop->getID()] = prop,
707 <  //it will cause memory leak if we already contain a propery which has the same name of prop
708 <  
709 <  if(result != properties.end()){
703 > double SimInfo::calcMaxCutoffRadius() {
704 >
705 >
706 >    std::set<AtomType*> atomTypes;
707 >    std::set<AtomType*>::iterator i;
708 >    std::vector<double> cutoffRadius;
709 >
710 >    //get the unique atom types
711 >    atomTypes = getUniqueAtomTypes();
712 >
713 >    //query the max cutoff radius among these atom types
714 >    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
715 >        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
716 >    }
717 >
718 >    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
719 > #ifdef IS_MPI
720 >    //pick the max cutoff radius among the processors
721 > #endif
722 >
723 >    return maxCutoffRadius;
724 > }
725 >
726 > void SimInfo::setupCutoff() {
727 >    double rcut_;  //cutoff radius
728 >    double rsw_; //switching radius
729      
730 <    delete (*result).second;
731 <    (*result).second = prop;
730 >    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
731 >        
732 >        if (!simParams_->haveRcut()){
733 >            sprintf(painCave.errMsg,
734 >                "SimCreator Warning: No value was set for the cutoffRadius.\n"
735 >                "\tOOPSE will use a default value of 15.0 angstroms"
736 >                "\tfor the cutoffRadius.\n");
737 >            painCave.isFatal = 0;
738 >            simError();
739 >            rcut_ = 15.0;
740 >        } else{
741 >            rcut_ = simParams_->getRcut();
742 >        }
743 >
744 >        if (!simParams_->haveRsw()){
745 >            sprintf(painCave.errMsg,
746 >                "SimCreator Warning: No value was set for switchingRadius.\n"
747 >                "\tOOPSE will use a default value of\n"
748 >                "\t0.95 * cutoffRadius for the switchingRadius\n");
749 >            painCave.isFatal = 0;
750 >            simError();
751 >            rsw_ = 0.95 * rcut_;
752 >        } else{
753 >            rsw_ = simParams_->getRsw();
754 >        }
755 >
756 >    } else {
757 >        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
758 >        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
759 >        
760 >        if (simParams_->haveRcut()) {
761 >            rcut_ = simParams_->getRcut();
762 >        } else {
763 >            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
764 >            rcut_ = calcMaxCutoffRadius();
765 >        }
766 >
767 >        if (simParams_->haveRsw()) {
768 >            rsw_  = simParams_->getRsw();
769 >        } else {
770 >            rsw_ = rcut_;
771 >        }
772 >    
773 >    }
774 >        
775 >    double rnblist = rcut_ + 1; // skin of neighbor list
776 >
777 >    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
778 >    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
779 > }
780 >
781 > void SimInfo::addProperty(GenericData* genData) {
782 >    properties_.addProperty(genData);  
783 > }
784 >
785 > void SimInfo::removeProperty(const std::string& propName) {
786 >    properties_.removeProperty(propName);  
787 > }
788 >
789 > void SimInfo::clearProperties() {
790 >    properties_.clearProperties();
791 > }
792 >
793 > std::vector<std::string> SimInfo::getPropertyNames() {
794 >    return properties_.getPropertyNames();  
795 > }
796        
797 <  }
798 <  else{
797 > std::vector<GenericData*> SimInfo::getProperties() {
798 >    return properties_.getProperties();
799 > }
800  
801 <    properties[prop->getID()] = prop;
801 > GenericData* SimInfo::getPropertyByName(const std::string& propName) {
802 >    return properties_.getPropertyByName(propName);
803 > }
804  
805 <  }
805 > void SimInfo::setSnapshotManager(SnapshotManager* sman) {
806 >    sman_ = sman;
807 >
808 >    Molecule* mol;
809 >    RigidBody* rb;
810 >    Atom* atom;
811 >    SimInfo::MoleculeIterator mi;
812 >    Molecule::RigidBodyIterator rbIter;
813 >    Molecule::AtomIterator atomIter;;
814 >
815 >    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
816 >        
817 >        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
818 >            atom->setSnapshotManager(sman_);
819 >        }
820 >        
821 >        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
822 >            rb->setSnapshotManager(sman_);
823 >        }
824 >    }    
825      
826   }
827  
828 < GenericData* SimInfo::getProperty(const string& propName){
828 > Vector3d SimInfo::getComVel(){
829 >    SimInfo::MoleculeIterator i;
830 >    Molecule* mol;
831 >
832 >    Vector3d comVel(0.0);
833 >    double totalMass = 0.0;
834 >    
835  
836 <  map<string, GenericData*>::iterator result;
837 <  
838 <  //string lowerCaseName = ();
839 <  
840 <  result = properties.find(propName);
841 <  
842 <  if(result != properties.end())
843 <    return (*result).second;  
844 <  else  
571 <    return NULL;  
836 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
837 >        double mass = mol->getMass();
838 >        totalMass += mass;
839 >        comVel += mass * mol->getComVel();
840 >    }  
841 >
842 >    comVel /= totalMass;
843 >
844 >    return comVel;
845   }
846  
847 + Vector3d SimInfo::getCom(){
848 +    SimInfo::MoleculeIterator i;
849 +    Molecule* mol;
850  
851 < void SimInfo::getFortranGroupArrays(SimInfo* info,
852 <                                    vector<int>& FglobalGroupMembership,
853 <                                    vector<double>& mfact){
854 <  
855 <  Molecule* myMols;
856 <  Atom** myAtoms;
857 <  int numAtom;
858 <  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 <  
851 >    Vector3d com(0.0);
852 >    double totalMass = 0.0;
853 >    
854 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
855 >        double mass = mol->getMass();
856 >        totalMass += mass;
857 >        com += mass * mol->getCom();
858 >    }  
859  
860 <  // Fix the silly fortran indexing problem
597 < #ifdef IS_MPI
598 <  numAtom = mpiSim->getNAtomsGlobal();
599 < #else
600 <  numAtom = n_atoms;
601 < #endif
602 <  for (int i = 0; i < numAtom; i++)
603 <    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
604 <  
860 >    com /= totalMass;
861  
862 <  myMols = info->molecules;
607 <  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)){
862 >    return com;
863  
864 <      totalMass = myCutoffGroup->getMass();
615 <      
616 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
617 <          cutoffAtom != NULL;
618 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
619 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
620 <      }  
621 <    }
622 <  }
864 > }        
865  
866 + std::ostream& operator <<(std::ostream& o, SimInfo& info) {
867 +
868 +    return o;
869   }
870 +
871 + }//end namespace oopse
872 +

Diff Legend

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