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

Comparing:
trunk/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1617 by chuckv, Wed Oct 20 20:46:20 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1727 by tim, Thu Nov 11 16:41:58 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 +
35   #include "brains/SimInfo.hpp"
36 < #define __C
10 < #include "brains/fSimulation.h"
11 < #include "utils/simError.h"
12 < #include "UseTheForce/DarkSide/simulation_interface.h"
13 < #include "UseTheForce/notifyCutoffs_interface.h"
36 > #include "utils/MemoryUtils.hpp"
37  
38 < //#include "UseTheForce/fortranWrappers.hpp"
38 > namespace oopse {
39  
40 < #include "math/MatVec3.h"
40 > SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
41 >        nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){
42  
19 – #ifdef IS_MPI
20 – #include "brains/mpiSimulation.hpp"
21 – #endif
22 –
23 – inline double roundMe( double x ){
24 –  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
43   }
26 –          
27 – inline double min( double a, double b ){
28 –  return (a < b ) ? a : b;
29 – }
44  
45 < SimInfo* currentInfo;
45 > SimInfo::~SimInfo() {
46 >    //MemoryUtils::deleteVectorOfPointer(molecules_);
47 >    delete sman_;
48  
49 < SimInfo::SimInfo(){
49 > }
50  
35 –  n_constraints = 0;
36 –  nZconstraints = 0;
37 –  n_oriented = 0;
38 –  n_dipoles = 0;
39 –  ndf = 0;
40 –  ndfRaw = 0;
41 –  nZconstraints = 0;
42 –  the_integrator = NULL;
43 –  setTemp = 0;
44 –  thermalTime = 0.0;
45 –  currentTime = 0.0;
46 –  rCut = 0.0;
47 –  rSw = 0.0;
51  
52 <  haveRcut = 0;
53 <  haveRsw = 0;
54 <  boxIsInit = 0;
55 <  
53 <  resetTime = 1e99;
52 > bool SimInfo::addMolecule(Molecule* mol) {
53 >    MoleculeIterator i;
54 >    i = std::find(molecules_.begin(), molecules_.end(), mol);
55 >    if (i != molecules_.end() ) {
56  
57 <  orthoRhombic = 0;
58 <  orthoTolerance = 1E-6;
59 <  useInitXSstate = true;
57 >        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
58 >        
59 >        nAtoms_ += mol->getNAtoms();
60 >        nBonds_ += mol->getNBonds();
61 >        nBends_ += mol->getNBends();
62 >        nTorsions_ += mol->getNTorsions();
63 >        nRigidBodies_ += mol->getNRigidBodies();
64 >        nIntegrableObjects_ += mol->getNIntegrableObjects();
65 >        nCutoffGroups_ += mol->getNCutoffGroups();
66 >        nConstraints_ += mol->getNConstraints();
67  
68 <  usePBC = 0;
69 <  useLJ = 0;
70 <  useSticky = 0;
71 <  useCharges = 0;
72 <  useDipoles = 0;
64 <  useReactionField = 0;
65 <  useGB = 0;
66 <  useEAM = 0;
67 <  useSolidThermInt = 0;
68 <  useLiquidThermInt = 0;
68 >        return true;
69 >    } else {
70 >        return false;
71 >    }
72 > }
73  
74 <  haveCutoffGroups = false;
74 > bool SimInfo::removeMolecule(Molecule* mol) {
75 >    MoleculeIterator i;
76 >    i = std::find(molecules_.begin(), molecules_.end(), mol);
77  
78 <  excludes = Exclude::Instance();
78 >    if (i != molecules_.end() ) {
79  
80 <  myConfiguration = new SimState();
80 >        nAtoms_ -= mol->getNAtoms();
81 >        nBonds_ -= mol->getNBonds();
82 >        nBends_ -= mol->getNBends();
83 >        nTorsions_ -= mol->getNTorsions();
84 >        nRigidBodies_ -= mol->getNRigidBodies();
85 >        nIntegrableObjects_ -= mol->getNIntegrableObjects();
86 >        nCutoffGroups_ -= mol->getNCutoffGroups();
87 >        nConstraints_ -= mol->getNConstraints();
88  
89 <  has_minimizer = false;
77 <  the_minimizer =NULL;
89 >        molecules_.erase(mol->getGlobalIndex());
90  
91 <  ngroup = 0;
91 >        delete mol;
92 >        
93 >        return true;
94 >    } else {
95 >        return false;
96 >    }
97  
81 – }
98  
99 + }    
100  
101 < SimInfo::~SimInfo(){
101 >        
102 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
103 >    i = molecules_.begin();
104 >    return i == molecules_.end() ? NULL : *i;
105 > }    
106  
107 <  delete myConfiguration;
107 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
108 >    ++i;
109 >    return i == molecules_.end() ? NULL : *i;    
110 > }
111  
88 –  map<string, GenericData*>::iterator i;
89 –  
90 –  for(i = properties.begin(); i != properties.end(); i++)
91 –    delete (*i).second;
112  
113 < }
113 > void SimInfo::calcNdf() {
114 >    int ndf_local;
115 >    MoleculeIterator i;
116 >    std::vector<StuntDouble*>::iterator j;
117 >    Molecule* mol;
118 >    StuntDouble* integrableObject;
119  
120 < void SimInfo::setBox(double newBox[3]) {
121 <  
122 <  int i, j;
123 <  double tempMat[3][3];
120 >    ndf_local = 0;
121 >    
122 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
123 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
124 >               integrableObject = mol->nextIntegrableObject(j)) {
125  
126 <  for(i=0; i<3; i++)
101 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
126 >            ndf_local += 3;
127  
128 <  tempMat[0][0] = newBox[0];
129 <  tempMat[1][1] = newBox[1];
130 <  tempMat[2][2] = newBox[2];
128 >            if (integrableObject->isDirectional()) {
129 >                if (integrableObject->isLinear()) {
130 >                    ndf_local += 2;
131 >                } else {
132 >                    ndf_local += 3;
133 >                }
134 >            }
135 >            
136 >        }//end for (integrableObject)
137 >    }// end for (mol)
138 >    
139 >    // n_constraints is local, so subtract them on each processor
140 >    ndf_local -= nConstraints_;
141  
142 <  setBoxM( tempMat );
142 > #ifdef IS_MPI
143 >    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
144 > #else
145 >    ndf_ = ndf_local;
146 > #endif
147  
148 +    // nZconstraints is global, as are the 3 COM translations for the
149 +    // entire system:
150 +    ndf_ = ndf_ - 3 - nZconstraints;
151 +
152   }
153  
154 < void SimInfo::setBoxM( double theBox[3][3] ){
155 <  
113 <  int i, j;
114 <  double FortranHmat[9]; // to preserve compatibility with Fortran the
115 <                         // ordering in the array is as follows:
116 <                         // [ 0 3 6 ]
117 <                         // [ 1 4 7 ]
118 <                         // [ 2 5 8 ]
119 <  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
154 > void SimInfo::calcNdfRaw() {
155 >    int ndfRaw_local;
156  
157 <  if( !boxIsInit ) boxIsInit = 1;
157 >    MoleculeIterator i;
158 >    std::vector<StuntDouble*>::iterator j;
159 >    Molecule* mol;
160 >    StuntDouble* integrableObject;
161  
162 <  for(i=0; i < 3; i++)
163 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
164 <  
165 <  calcBoxL();
166 <  calcHmatInv();
162 >    // Raw degrees of freedom that we have to set
163 >    ndfRaw_local = 0;
164 >    
165 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
166 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
167 >               integrableObject = mol->nextIntegrableObject(j)) {
168  
169 <  for(i=0; i < 3; i++) {
130 <    for (j=0; j < 3; j++) {
131 <      FortranHmat[3*j + i] = Hmat[i][j];
132 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
133 <    }
134 <  }
169 >            ndfRaw_local += 3;
170  
171 <  setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic);
172 <
171 >            if (integrableObject->isDirectional()) {
172 >                if (integrableObject->isLinear()) {
173 >                    ndfRaw_local += 2;
174 >                } else {
175 >                    ndfRaw_local += 3;
176 >                }
177 >            }
178 >            
179 >        }
180 >    }
181 >    
182 > #ifdef IS_MPI
183 >    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
184 > #else
185 >    ndfRaw_ = ndfRaw_local;
186 > #endif
187   }
139 –
188  
189 < void SimInfo::getBoxM (double theBox[3][3]) {
189 > void SimInfo::calcNdfTrans() {
190 >    int ndfTrans_local;
191  
192 <  int i, j;
144 <  for(i=0; i<3; i++)
145 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
146 < }
192 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
193  
194  
195 < void SimInfo::scaleBox(double scale) {
196 <  double theBox[3][3];
197 <  int i, j;
195 > #ifdef IS_MPI
196 >    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
197 > #else
198 >    ndfTrans_ = ndfTrans_local;
199 > #endif
200  
201 <  // cerr << "Scaling box by " << scale << "\n";
202 <
155 <  for(i=0; i<3; i++)
156 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
157 <
158 <  setBoxM(theBox);
159 <
201 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
202 >
203   }
204  
205 < void SimInfo::calcHmatInv( void ) {
206 <  
207 <  int oldOrtho;
208 <  int i,j;
209 <  double smallDiag;
210 <  double tol;
211 <  double sanity[3][3];
205 > void SimInfo::addExcludePairs(Molecule* mol) {
206 >    std::vector<Bond*>::iterator bondIter;
207 >    std::vector<Bend*>::iterator bendIter;
208 >    std::vector<Torsion*>::iterator torsionIter;
209 >    Bond* bond;
210 >    Bend* bend;
211 >    Torsion* torsion;
212 >    int a;
213 >    int b;
214 >    int c;
215 >    int d;
216 >    
217 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
218 >        a = bond->getAtomA()->getGlobalIndex();
219 >        b = bond->getAtomB()->getGlobalIndex();        
220 >        exclude_.addPair(a, b);
221 >    }
222  
223 <  invertMat3( Hmat, HmatInv );
223 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
224 >        a = bend->getAtomA()->getGlobalIndex();
225 >        b = bend->getAtomB()->getGlobalIndex();        
226 >        c = bend->getAtomC()->getGlobalIndex();
227  
228 <  // check to see if Hmat is orthorhombic
229 <  
230 <  oldOrtho = orthoRhombic;
228 >        exclude_.addPair(a, b);
229 >        exclude_.addPair(a, c);
230 >        exclude_.addPair(b, c);        
231 >    }
232  
233 <  smallDiag = fabs(Hmat[0][0]);
234 <  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
235 <  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
236 <  tol = smallDiag * orthoTolerance;
233 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
234 >        a = torsion->getAtomA()->getGlobalIndex();
235 >        b = torsion->getAtomB()->getGlobalIndex();        
236 >        c = torsion->getAtomC()->getGlobalIndex();        
237 >        d = torsion->getAtomD()->getGlobalIndex();        
238  
239 <  orthoRhombic = 1;
240 <  
241 <  for (i = 0; i < 3; i++ ) {
242 <    for (j = 0 ; j < 3; j++) {
243 <      if (i != j) {
244 <        if (orthoRhombic) {
187 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
188 <        }        
189 <      }
239 >        exclude_.addPair(a, b);
240 >        exclude_.addPair(a, c);
241 >        exclude_.addPair(a, d);
242 >        exclude_.addPair(b, c);
243 >        exclude_.addPair(b, d);
244 >        exclude_.addPair(c, d);        
245      }
191 –  }
246  
193 –  if( oldOrtho != orthoRhombic ){
247      
195 –    if( orthoRhombic ) {
196 –      sprintf( painCave.errMsg,
197 –               "OOPSE is switching from the default Non-Orthorhombic\n"
198 –               "\tto the faster Orthorhombic periodic boundary computations.\n"
199 –               "\tThis is usually a good thing, but if you wan't the\n"
200 –               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
201 –               "\tvariable ( currently set to %G ) smaller.\n",
202 –               orthoTolerance);
203 –      painCave.severity = OOPSE_INFO;
204 –      simError();
205 –    }
206 –    else {
207 –      sprintf( painCave.errMsg,
208 –               "OOPSE is switching from the faster Orthorhombic to the more\n"
209 –               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
210 –               "\tThis is usually because the box has deformed under\n"
211 –               "\tNPTf integration. If you wan't to live on the edge with\n"
212 –               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
213 –               "\tvariable ( currently set to %G ) larger.\n",
214 –               orthoTolerance);
215 –      painCave.severity = OOPSE_WARNING;
216 –      simError();
217 –    }
218 –  }
248   }
249  
250 < void SimInfo::calcBoxL( void ){
251 <
252 <  double dx, dy, dz, dsq;
253 <
254 <  // boxVol = Determinant of Hmat
255 <
256 <  boxVol = matDet3( Hmat );
257 <
258 <  // boxLx
259 <  
260 <  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
232 <  dsq = dx*dx + dy*dy + dz*dz;
233 <  boxL[0] = sqrt( dsq );
234 <  //maxCutoff = 0.5 * boxL[0];
235 <
236 <  // boxLy
237 <  
238 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
239 <  dsq = dx*dx + dy*dy + dz*dz;
240 <  boxL[1] = sqrt( dsq );
241 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
242 <
243 <
244 <  // boxLz
245 <  
246 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
247 <  dsq = dx*dx + dy*dy + dz*dz;
248 <  boxL[2] = sqrt( dsq );
249 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
250 <
251 <  //calculate the max cutoff
252 <  maxCutoff =  calcMaxCutOff();
253 <  
254 <  checkCutOffs();
255 <
256 < }
257 <
258 <
259 < double SimInfo::calcMaxCutOff(){
260 <
261 <  double ri[3], rj[3], rk[3];
262 <  double rij[3], rjk[3], rki[3];
263 <  double minDist;
264 <
265 <  ri[0] = Hmat[0][0];
266 <  ri[1] = Hmat[1][0];
267 <  ri[2] = Hmat[2][0];
268 <
269 <  rj[0] = Hmat[0][1];
270 <  rj[1] = Hmat[1][1];
271 <  rj[2] = Hmat[2][1];
272 <
273 <  rk[0] = Hmat[0][2];
274 <  rk[1] = Hmat[1][2];
275 <  rk[2] = Hmat[2][2];
276 <    
277 <  crossProduct3(ri, rj, rij);
278 <  distXY = dotProduct3(rk,rij) / norm3(rij);
279 <
280 <  crossProduct3(rj,rk, rjk);
281 <  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
282 <
283 <  crossProduct3(rk,ri, rki);
284 <  distZX = dotProduct3(rj,rki) / norm3(rki);
285 <
286 <  minDist = min(min(distXY, distYZ), distZX);
287 <  return minDist/2;
288 <  
289 < }
290 <
291 < void SimInfo::wrapVector( double thePos[3] ){
292 <
293 <  int i;
294 <  double scaled[3];
295 <
296 <  if( !orthoRhombic ){
297 <    // calc the scaled coordinates.
298 <  
299 <
300 <    matVecMul3(HmatInv, thePos, scaled);
250 > void SimInfo::removeExcludePairs(Molecule* mol) {
251 >    std::vector<Bond*>::iterator bondIter;
252 >    std::vector<Bend*>::iterator bendIter;
253 >    std::vector<Torsion*>::iterator torsionIter;
254 >    Bond* bond;
255 >    Bend* bend;
256 >    Torsion* torsion;
257 >    int a;
258 >    int b;
259 >    int c;
260 >    int d;
261      
262 <    for(i=0; i<3; i++)
263 <      scaled[i] -= roundMe(scaled[i]);
264 <    
265 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
266 <    
307 <    matVecMul3(Hmat, scaled, thePos);
262 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
263 >        a = bond->getAtomA()->getGlobalIndex();
264 >        b = bond->getAtomB()->getGlobalIndex();        
265 >        exclude_.removePair(a, b);
266 >    }
267  
268 <  }
269 <  else{
270 <    // calc the scaled coordinates.
271 <    
313 <    for(i=0; i<3; i++)
314 <      scaled[i] = thePos[i]*HmatInv[i][i];
315 <    
316 <    // wrap the scaled coordinates
317 <    
318 <    for(i=0; i<3; i++)
319 <      scaled[i] -= roundMe(scaled[i]);
320 <    
321 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
322 <    
323 <    for(i=0; i<3; i++)
324 <      thePos[i] = scaled[i]*Hmat[i][i];
325 <  }
326 <    
327 < }
268 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
269 >        a = bend->getAtomA()->getGlobalIndex();
270 >        b = bend->getAtomB()->getGlobalIndex();        
271 >        c = bend->getAtomC()->getGlobalIndex();
272  
273 <
274 < int SimInfo::getNDF(){
275 <  int ndf_local;
332 <
333 <  ndf_local = 0;
334 <  
335 <  for(int i = 0; i < integrableObjects.size(); i++){
336 <    ndf_local += 3;
337 <    if (integrableObjects[i]->isDirectional()) {
338 <      if (integrableObjects[i]->isLinear())
339 <        ndf_local += 2;
340 <      else
341 <        ndf_local += 3;
273 >        exclude_.removePair(a, b);
274 >        exclude_.removePair(a, c);
275 >        exclude_.removePair(b, c);        
276      }
343 –  }
277  
278 <  // n_constraints is local, so subtract them on each processor:
278 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
279 >        a = torsion->getAtomA()->getGlobalIndex();
280 >        b = torsion->getAtomB()->getGlobalIndex();        
281 >        c = torsion->getAtomC()->getGlobalIndex();        
282 >        d = torsion->getAtomD()->getGlobalIndex();        
283  
284 <  ndf_local -= n_constraints;
285 <
286 < #ifdef IS_MPI
287 <  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
288 < #else
289 <  ndf = ndf_local;
353 < #endif
354 <
355 <  // nZconstraints is global, as are the 3 COM translations for the
356 <  // entire system:
357 <
358 <  ndf = ndf - 3 - nZconstraints;
359 <
360 <  return ndf;
361 < }
362 <
363 < int SimInfo::getNDFraw() {
364 <  int ndfRaw_local;
365 <
366 <  // Raw degrees of freedom that we have to set
367 <  ndfRaw_local = 0;
368 <
369 <  for(int i = 0; i < integrableObjects.size(); i++){
370 <    ndfRaw_local += 3;
371 <    if (integrableObjects[i]->isDirectional()) {
372 <       if (integrableObjects[i]->isLinear())
373 <        ndfRaw_local += 2;
374 <      else
375 <        ndfRaw_local += 3;
284 >        exclude_.removePair(a, b);
285 >        exclude_.removePair(a, c);
286 >        exclude_.removePair(a, d);
287 >        exclude_.removePair(b, c);
288 >        exclude_.removePair(b, d);
289 >        exclude_.removePair(c, d);        
290      }
377 –  }
378 –    
379 – #ifdef IS_MPI
380 –  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
381 – #else
382 –  ndfRaw = ndfRaw_local;
383 – #endif
291  
385 –  return ndfRaw;
292   }
293  
388 – int SimInfo::getNDFtranslational() {
389 –  int ndfTrans_local;
294  
295 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
295 > void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
296 >    int curStampId;
297  
298 +    //index from 0
299 +    curStampId = molStampIds_.size();
300  
301 < #ifdef IS_MPI
302 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
396 < #else
397 <  ndfTrans = ndfTrans_local;
398 < #endif
399 <
400 <  ndfTrans = ndfTrans - 3 - nZconstraints;
401 <
402 <  return ndfTrans;
301 >    moleculeStamps_.push_back(molStamp);
302 >    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
303   }
304  
305 < int SimInfo::getTotIntegrableObjects() {
406 <  int nObjs_local;
407 <  int nObjs;
305 > void SimInfo::update() {
306  
409 –  nObjs_local =  integrableObjects.size();
307  
411 –
412 – #ifdef IS_MPI
413 –  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
414 – #else
415 –  nObjs = nObjs_local;
416 – #endif
417 –
418 –
419 –  return nObjs;
420 – }
421 –
422 – void SimInfo::refreshSim(){
423 –
424 –  simtype fInfo;
425 –  int isError;
426 –  int n_global;
427 –  int* excl;
428 –
429 –  fInfo.dielect = 0.0;
430 –
431 –  if( useDipoles ){
432 –    if( useReactionField )fInfo.dielect = dielectric;
433 –  }
434 –
435 –  fInfo.SIM_uses_PBC = usePBC;
436 –  //fInfo.SIM_uses_LJ = 0;
437 –  fInfo.SIM_uses_LJ = useLJ;
438 –  fInfo.SIM_uses_sticky = useSticky;
439 –  //fInfo.SIM_uses_sticky = 0;
440 –  fInfo.SIM_uses_charges = useCharges;
441 –  fInfo.SIM_uses_dipoles = useDipoles;
442 –  //fInfo.SIM_uses_dipoles = 0;
443 –  fInfo.SIM_uses_RF = useReactionField;
444 –  //fInfo.SIM_uses_RF = 0;
445 –  fInfo.SIM_uses_GB = useGB;
446 –  fInfo.SIM_uses_EAM = useEAM;
447 –
448 –  n_exclude = excludes->getSize();
449 –  excl = excludes->getFortranArray();
450 –  
451 – #ifdef IS_MPI
452 –  n_global = mpiSim->getNAtomsGlobal();
453 – #else
454 –  n_global = n_atoms;
455 – #endif
456 –  
457 –  isError = 0;
458 –  
459 –  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
460 –  //it may not be a good idea to pass the address of first element in vector
461 –  //since c++ standard does not require vector to be stored continuously in meomory
462 –  //Most of the compilers will organize the memory of vector continuously
463 –  setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
464 –                  &nGlobalExcludes, globalExcludes, molMembershipArray,
465 –                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
466 –
467 –  if( isError ){
308      
309 <    sprintf( painCave.errMsg,
310 <             "There was an error setting the simulation information in fortran.\n" );
311 <    painCave.isFatal = 1;
312 <    painCave.severity = OOPSE_ERROR;
313 <    simError();
314 <  }
315 <  
316 < #ifdef IS_MPI
317 <  sprintf( checkPointMsg,
318 <           "succesfully sent the simulation information to fortran.\n");
319 <  MPIcheckPoint();
480 < #endif // is_mpi
481 <  
482 <  this->ndf = this->getNDF();
483 <  this->ndfRaw = this->getNDFraw();
484 <  this->ndfTrans = this->getNDFtranslational();
485 < }
309 >    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
310 >    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
311 >    std::vector<int> localToGlobalCutoffGroupIndex;
312 >    typename SimInfo::MoleculeIterator mi;
313 >    typename Molecule::AtomIterator ai;
314 >    typename Molecule::CutoffGroupIterator ci;
315 >    Molecule* mol;
316 >    Atom* atom;
317 >    CutoffGroup* cg;
318 >    mpiSimData parallelData;
319 >    int isError;
320  
321 < void SimInfo::setDefaultRcut( double theRcut ){
488 <  
489 <  haveRcut = 1;
490 <  rCut = theRcut;
491 <  rList = rCut + 1.0;
492 <  
493 <  notifyFortranCutoffs( &rCut, &rSw, &rList );
494 < }
321 >    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
322  
323 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
323 >        //local index(index in DataStorge) of atom is important
324 >        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
325 >            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
326 >        }
327  
328 <  rSw = theRsw;
329 <  setDefaultRcut( theRcut );
330 < }
328 >        //local index of cutoff group is trivial, it only depends on the order of travesing
329 >        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
330 >            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
331 >        }        
332 >        
333 >    }
334  
335 <
336 < void SimInfo::checkCutOffs( void ){
337 <  
338 <  if( boxIsInit ){
335 >    //Setup Parallel Data and pass the index arrays to fortran
336 >    parallelData.nMolGlobal = getNMolGlobal();
337 >    parallelData.nMolLocal = ;
338 >    parallelData.nAtomsGlobal = ;
339 >    parallelData.nAtomsLocal = ;
340 >    parallelData.nGroupsGlobal = ;
341 >    parallelData.nGroupsLocal = ;
342 >    parallelData.myNode = worldRank;
343 >    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
344      
345 <    //we need to check cutOffs against the box
346 <    
347 <    if( rCut > maxCutoff ){
510 <      sprintf( painCave.errMsg,
511 <               "cutoffRadius is too large for the current periodic box.\n"
512 <               "\tCurrent Value of cutoffRadius = %G at time %G\n "
513 <               "\tThis is larger than half of at least one of the\n"
514 <               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
515 <               "\n"
516 <               "\t[ %G %G %G ]\n"
517 <               "\t[ %G %G %G ]\n"
518 <               "\t[ %G %G %G ]\n",
519 <               rCut, currentTime,
520 <               Hmat[0][0], Hmat[0][1], Hmat[0][2],
521 <               Hmat[1][0], Hmat[1][1], Hmat[1][2],
522 <               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
523 <      painCave.severity = OOPSE_ERROR;
524 <      painCave.isFatal = 1;
525 <      simError();
526 <    }    
527 <  } else {
528 <    // initialize this stuff before using it, OK?
529 <    sprintf( painCave.errMsg,
530 <             "Trying to check cutoffs without a box.\n"
531 <             "\tOOPSE should have better programmers than that.\n" );
532 <    painCave.severity = OOPSE_ERROR;
533 <    painCave.isFatal = 1;
534 <    simError();      
535 <  }
536 <  
537 < }
345 >    setFsimParallel(parallelData,            &(parallelData->nAtomsLocal),
346 >                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
347 >                    &localToGlobalCutoffGroupIndex[0], &isError);
348  
349 < void SimInfo::addProperty(GenericData* prop){
349 >    if (isError) {
350 >        sprintf(painCave.errMsg,
351 >                "mpiRefresh errror: fortran didn't like something we gave it.\n");
352 >        painCave.isFatal = 1;
353 >        simError();
354 >    }
355  
356 <  map<string, GenericData*>::iterator result;
357 <  result = properties.find(prop->getID());
543 <  
544 <  //we can't simply use  properties[prop->getID()] = prop,
545 <  //it will cause memory leak if we already contain a propery which has the same name of prop
546 <  
547 <  if(result != properties.end()){
548 <    
549 <    delete (*result).second;
550 <    (*result).second = prop;
551 <      
552 <  }
553 <  else{
356 >    sprintf(checkPointMsg, " mpiRefresh successful.\n");
357 >    MPIcheckPoint();
358  
555 –    properties[prop->getID()] = prop;
359  
557 –  }
558 –    
360   }
361  
362 < GenericData* SimInfo::getProperty(const string& propName){
562 <
563 <  map<string, GenericData*>::iterator result;
564 <  
565 <  //string lowerCaseName = ();
566 <  
567 <  result = properties.find(propName);
568 <  
569 <  if(result != properties.end())
570 <    return (*result).second;  
571 <  else  
572 <    return NULL;  
573 < }
362 > std::ostream& operator <<(ostream& o, SimInfo& info) {
363  
364 <
576 < void SimInfo::getFortranGroupArrays(SimInfo* info,
577 <                                    vector<int>& FglobalGroupMembership,
578 <                                    vector<double>& mfact){
579 <  
580 <  Molecule* myMols;
581 <  Atom** myAtoms;
582 <  int numAtom;
583 <  double mtot;
584 <  int numMol;
585 <  int numCutoffGroups;
586 <  CutoffGroup* myCutoffGroup;
587 <  vector<CutoffGroup*>::iterator iterCutoff;
588 <  Atom* cutoffAtom;
589 <  vector<Atom*>::iterator iterAtom;
590 <  int atomIndex;
591 <  double totalMass;
592 <  
593 <  mfact.clear();
594 <  FglobalGroupMembership.clear();
595 <  
596 <
597 <  // Fix the silly fortran indexing problem
598 < #ifdef IS_MPI
599 <  numAtom = mpiSim->getNAtomsGlobal();
600 < #else
601 <  numAtom = n_atoms;
602 < #endif
603 <  for (int i = 0; i < numAtom; i++)
604 <    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
605 <  
606 <
607 <  myMols = info->molecules;
608 <  numMol = info->n_mol;
609 <  for(int i  = 0; i < numMol; i++){
610 <    numCutoffGroups = myMols[i].getNCutoffGroups();
611 <    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
612 <        myCutoffGroup != NULL;
613 <        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
614 <
615 <      totalMass = myCutoffGroup->getMass();
616 <      
617 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
618 <          cutoffAtom != NULL;
619 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
620 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
621 <      }  
622 <    }
623 <  }
624 <
364 >    return o;
365   }
366 +
367 + }//end namespace oopse

Diff Legend

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