ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/brains/SimInfo.cpp
(Generate patch)

Comparing trunk/src/brains/SimInfo.cpp (file contents):
Revision 143 by chrisfen, Fri Oct 22 22:54:01 2004 UTC vs.
Revision 328 by tim, Sun Feb 13 20:36:24 2005 UTC

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

Diff Legend

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