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

Diff Legend

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