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

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

Diff Legend

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