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

Comparing:
trunk/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1490 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents), Revision 1710 by tim, Thu Nov 4 19:48:22 2004 UTC

# Line 1 | Line 1
1 < #include <stdlib.h>
2 < #include <string.h>
3 < #include <math.h>
1 > /*
2 > * Copyright (C) 2000-2004  Object Oriented Parallel Simulation Engine (OOPSE) project
3 > *
4 > * Contact: oopse@oopse.org
5 > *
6 > * This program is free software; you can redistribute it and/or
7 > * modify it under the terms of the GNU Lesser General Public License
8 > * as published by the Free Software Foundation; either version 2.1
9 > * of the License, or (at your option) any later version.
10 > * All we ask is that proper credit is given for our work, which includes
11 > * - but is not limited to - adding the above copyright notice to the beginning
12 > * of your source code files, and to any copyright notice that you may distribute
13 > * with programs based on this work.
14 > *
15 > * This program is distributed in the hope that it will be useful,
16 > * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 > * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
18 > * GNU Lesser General Public License for more details.
19 > *
20 > * You should have received a copy of the GNU Lesser General Public License
21 > * along with this program; if not, write to the Free Software
22 > * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
23 > *
24 > */
25  
26 < #include <iostream>
27 < using namespace std;
26 > /**
27 > * @file SimInfo.cpp
28 > * @author    tlin
29 > * @date  11/02/2004
30 > * @version 1.0
31 > */
32  
33 < #include "SimInfo.hpp"
9 < #define __C
10 < #include "fSimulation.h"
11 < #include "simError.h"
33 > #include <algorithm>
34  
35 < #include "fortranWrappers.hpp"
35 > #include "brains/SimInfo.hpp"
36 > #include "utils/MemoryUtils.hpp"
37  
38 < #include "MatVec3.h"
38 > namespace oopse {
39  
40 < #ifdef IS_MPI
41 < #include "mpiSimulation.hpp"
19 < #endif
40 > SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
41 >        nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){
42  
21 – inline double roundMe( double x ){
22 –  return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 );
43   }
24 –          
25 – inline double min( double a, double b ){
26 –  return (a < b ) ? a : b;
27 – }
44  
45 < SimInfo* currentInfo;
45 > SimInfo::~SimInfo() {
46 >    MemoryUtils::deleteVectorOfPointer(molecules_);
47 >    delete sman_;
48  
31 – SimInfo::SimInfo(){
32 –
33 –  n_constraints = 0;
34 –  nZconstraints = 0;
35 –  n_oriented = 0;
36 –  n_dipoles = 0;
37 –  ndf = 0;
38 –  ndfRaw = 0;
39 –  nZconstraints = 0;
40 –  the_integrator = NULL;
41 –  setTemp = 0;
42 –  thermalTime = 0.0;
43 –  currentTime = 0.0;
44 –  rCut = 0.0;
45 –  rSw = 0.0;
46 –
47 –  haveRcut = 0;
48 –  haveRsw = 0;
49 –  boxIsInit = 0;
50 –  
51 –  resetTime = 1e99;
52 –
53 –  orthoRhombic = 0;
54 –  orthoTolerance = 1E-6;
55 –  useInitXSstate = true;
56 –
57 –  usePBC = 0;
58 –  useLJ = 0;
59 –  useSticky = 0;
60 –  useCharges = 0;
61 –  useDipoles = 0;
62 –  useReactionField = 0;
63 –  useGB = 0;
64 –  useEAM = 0;
65 –  useSolidThermInt = 0;
66 –  useLiquidThermInt = 0;
67 –
68 –  haveCutoffGroups = false;
69 –
70 –  excludes = Exclude::Instance();
71 –
72 –  myConfiguration = new SimState();
73 –
74 –  has_minimizer = false;
75 –  the_minimizer =NULL;
76 –
77 –  ngroup = 0;
78 –
79 –  wrapMeSimInfo( this );
49   }
50  
51  
52 < SimInfo::~SimInfo(){
52 > bool SimInfo::addMolecule(Molecule* mol) {
53 >    std::vector<Molecule*>::iterator i;
54 >    i = std::find(molecules_.begin(), molecules_.end(), mol);
55 >    if (i != molecules_.end() ) {
56 >        molecules_.push_back(mol);
57  
58 <  delete myConfiguration;
58 >        nAtoms_ += mol->getNAtoms();
59 >        nBonds_ += mol->getNBonds();
60 >        nBends_ += mol->getNBends();
61 >        nTorsions_ += mol->getNTorsions();
62 >        nRigidBodies_ += mol->getNRigidBodies();
63 >        nIntegrableObjects_ += mol->getNIntegrableObjects();
64 >        nCutoffGroups_ += mol->getNCutoffGroups();
65 >        nConstraints_ += mol->getNConstraints();
66  
67 <  map<string, GenericData*>::iterator i;
68 <  
69 <  for(i = properties.begin(); i != properties.end(); i++)
70 <    delete (*i).second;
91 <
67 >        return true;
68 >    } else {
69 >        return false;
70 >    }
71   }
72  
73 < void SimInfo::setBox(double newBox[3]) {
74 <  
75 <  int i, j;
97 <  double tempMat[3][3];
73 > bool SimInfo::removeMolecule(Molecule* mol) {
74 >    std::vector<Molecule*>::iterator i;
75 >    i = std::find(molecules_.begin(), molecules_.end(), mol);
76  
77 <  for(i=0; i<3; i++)
78 <    for (j=0; j<3; j++) tempMat[i][j] = 0.0;;
77 >    if (i != molecules_.end() ) {
78 >        molecules_.push_back(mol);
79 >        nAtoms_ -= mol->getNAtoms();
80 >        nBonds_ -= mol->getNBonds();
81 >        nBends_ -= mol->getNBends();
82 >        nTorsions_ -= mol->getNTorsions();
83 >        nRigidBodies_ -= mol->getNRigidBodies();
84 >        nIntegrableObjects_ -= mol->getNIntegrableObjects();
85 >        nCutoffGroups_ -= mol->getNCutoffGroups();
86 >        nConstraints_ -= mol->getNConstraints();
87  
88 <  tempMat[0][0] = newBox[0];
89 <  tempMat[1][1] = newBox[1];
90 <  tempMat[2][2] = newBox[2];
88 >        return true;
89 >    } else {
90 >        return false;
91 >    }
92  
106 –  setBoxM( tempMat );
93  
94 < }
94 > }    
95  
96 < void SimInfo::setBoxM( double theBox[3][3] ){
97 <  
98 <  int i, j;
99 <  double FortranHmat[9]; // to preserve compatibility with Fortran the
100 <                         // 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);
96 >        
97 > Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
98 >    i = molecules_.begin();
99 >    return i == molecules_.end() ? NULL : *i;
100 > }    
101  
102 <  if( !boxIsInit ) boxIsInit = 1;
103 <
104 <  for(i=0; i < 3; i++)
123 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
124 <  
125 <  calcBoxL();
126 <  calcHmatInv();
127 <
128 <  for(i=0; i < 3; i++) {
129 <    for (j=0; j < 3; j++) {
130 <      FortranHmat[3*j + i] = Hmat[i][j];
131 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
132 <    }
133 <  }
134 <
135 <  setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic);
136 <
102 > Molecule* SimInfo::nextMolecule(std::vector<Molecule*>::iterator& i) {
103 >    ++i;
104 >    return i == molecules_.end() ? NULL : *i;    
105   }
138 –
106  
140 – void SimInfo::getBoxM (double theBox[3][3]) {
107  
108 <  int i, j;
109 <  for(i=0; i<3; i++)
110 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
111 < }
108 > void SimInfo::calcNDF() {
109 >    int ndf_local;
110 >    std::vector<Molecule*>::iterator i;
111 >    std::vector<StuntDouble*>::iterator j;
112 >    Molecule* mol;
113 >    StuntDouble* integrableObject;
114  
115 +    ndf_local = 0;
116 +    
117 +    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
118 +        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
119 +               integrableObject = mol->nextIntegrableObject(j)) {
120  
121 < void SimInfo::scaleBox(double scale) {
149 <  double theBox[3][3];
150 <  int i, j;
121 >            ndf_local += 3;
122  
123 <  // cerr << "Scaling box by " << scale << "\n";
123 >            if (integrableObject->isDirectional()) {
124 >                if (integrableObject->isLinear()) {
125 >                    ndf_local += 2;
126 >                } else {
127 >                    ndf_local += 3;
128 >                }
129 >            }
130 >            
131 >        }//end for (integrableObject)
132 >    }// end for (mol)
133 >    
134 >    // n_constraints is local, so subtract them on each processor
135 >    ndf_local -= nConstraints_;
136  
137 <  for(i=0; i<3; i++)
138 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
137 > #ifdef IS_MPI
138 >    MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
139 > #else
140 >    ndf_ = ndf_local;
141 > #endif
142  
143 <  setBoxM(theBox);
143 >    // nZconstraints is global, as are the 3 COM translations for the
144 >    // entire system:
145 >    ndf_ = ndf_ - 3 - nZconstraints;
146  
147   }
148  
149 < void SimInfo::calcHmatInv( void ) {
150 <  
163 <  int oldOrtho;
164 <  int i,j;
165 <  double smallDiag;
166 <  double tol;
167 <  double sanity[3][3];
149 > void SimInfo::calcNDFRaw() {
150 >    int ndfRaw_local;
151  
152 <  invertMat3( Hmat, HmatInv );
152 >    std::vector<Molecule*>::iterator i;
153 >    std::vector<StuntDouble*>::iterator j;
154 >    Molecule* mol;
155 >    StuntDouble* integrableObject;
156  
157 <  // check to see if Hmat is orthorhombic
158 <  
159 <  oldOrtho = orthoRhombic;
157 >    // Raw degrees of freedom that we have to set
158 >    ndfRaw_local = 0;
159 >    
160 >    for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
161 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
162 >               integrableObject = mol->nextIntegrableObject(j)) {
163  
164 <  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;
164 >            ndfRaw_local += 3;
165  
166 <  orthoRhombic = 1;
167 <  
168 <  for (i = 0; i < 3; i++ ) {
169 <    for (j = 0 ; j < 3; j++) {
170 <      if (i != j) {
171 <        if (orthoRhombic) {
172 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
173 <        }        
174 <      }
166 >            if (integrableObject->isDirectional()) {
167 >                if (integrableObject->isLinear()) {
168 >                    ndfRaw_local += 2;
169 >                } else {
170 >                    ndfRaw_local += 3;
171 >                }
172 >            }
173 >            
174 >        }
175      }
190 –  }
191 –
192 –  if( oldOrtho != orthoRhombic ){
176      
177 <    if( orthoRhombic ) {
178 <      sprintf( painCave.errMsg,
179 <               "OOPSE is switching from the default Non-Orthorhombic\n"
180 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
181 <               "\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 <  }
177 > #ifdef IS_MPI
178 >    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
179 > #else
180 >    ndfRaw_ = ndfRaw_local;
181 > #endif
182   }
183  
184 < void SimInfo::calcBoxL( void ){
184 > void SimInfo::calcNDFTrans() {
185 >    int ndfTrans_local;
186  
187 <  double dx, dy, dz, dsq;
187 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
188  
224 –  // boxVol = Determinant of Hmat
189  
226 –  boxVol = matDet3( Hmat );
227 –
228 –  // boxLx
229 –  
230 –  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];
234 –
235 –  // boxLy
236 –  
237 –  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
238 –  dsq = dx*dx + dy*dy + dz*dz;
239 –  boxL[1] = sqrt( dsq );
240 –  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
241 –
242 –
243 –  // boxLz
244 –  
245 –  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
246 –  dsq = dx*dx + dy*dy + dz*dz;
247 –  boxL[2] = sqrt( dsq );
248 –  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
249 –
250 –  //calculate the max cutoff
251 –  maxCutoff =  calcMaxCutOff();
252 –  
253 –  checkCutOffs();
254 –
255 – }
256 –
257 –
258 – double SimInfo::calcMaxCutOff(){
259 –
260 –  double ri[3], rj[3], rk[3];
261 –  double rij[3], rjk[3], rki[3];
262 –  double minDist;
263 –
264 –  ri[0] = Hmat[0][0];
265 –  ri[1] = Hmat[1][0];
266 –  ri[2] = Hmat[2][0];
267 –
268 –  rj[0] = Hmat[0][1];
269 –  rj[1] = Hmat[1][1];
270 –  rj[2] = Hmat[2][1];
271 –
272 –  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);
278 –
279 –  crossProduct3(rj,rk, rjk);
280 –  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
281 –
282 –  crossProduct3(rk,ri, rki);
283 –  distZX = dotProduct3(rj,rki) / norm3(rki);
284 –
285 –  minDist = min(min(distXY, distYZ), distZX);
286 –  return minDist/2;
287 –  
288 – }
289 –
290 – void SimInfo::wrapVector( double thePos[3] ){
291 –
292 –  int i;
293 –  double scaled[3];
294 –
295 –  if( !orthoRhombic ){
296 –    // calc the scaled coordinates.
297 –  
298 –
299 –    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);
307 –
308 –  }
309 –  else{
310 –    // calc the scaled coordinates.
311 –    
312 –    for(i=0; i<3; i++)
313 –      scaled[i] = thePos[i]*HmatInv[i][i];
314 –    
315 –    // wrap the scaled coordinates
316 –    
317 –    for(i=0; i<3; i++)
318 –      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 –    
326 – }
327 –
328 –
329 – int SimInfo::getNDF(){
330 –  int ndf_local;
331 –
332 –  ndf_local = 0;
333 –  
334 –  for(int i = 0; i < integrableObjects.size(); i++){
335 –    ndf_local += 3;
336 –    if (integrableObjects[i]->isDirectional()) {
337 –      if (integrableObjects[i]->isLinear())
338 –        ndf_local += 2;
339 –      else
340 –        ndf_local += 3;
341 –    }
342 –  }
343 –
344 –  // n_constraints is local, so subtract them on each processor:
345 –
346 –  ndf_local -= n_constraints;
347 –
190   #ifdef IS_MPI
191 <  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
191 >    MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
192   #else
193 <  ndf = ndf_local;
193 >    ndfTrans_ = ndfTrans_local;
194   #endif
195  
196 <  // nZconstraints is global, as are the 3 COM translations for the
355 <  // entire system:
356 <
357 <  ndf = ndf - 3 - nZconstraints;
358 <
359 <  return ndf;
360 < }
361 <
362 < int SimInfo::getNDFraw() {
363 <  int ndfRaw_local;
364 <
365 <  // Raw degrees of freedom that we have to set
366 <  ndfRaw_local = 0;
367 <
368 <  for(int i = 0; i < integrableObjects.size(); i++){
369 <    ndfRaw_local += 3;
370 <    if (integrableObjects[i]->isDirectional()) {
371 <       if (integrableObjects[i]->isLinear())
372 <        ndfRaw_local += 2;
373 <      else
374 <        ndfRaw_local += 3;
375 <    }
376 <  }
377 <    
378 < #ifdef IS_MPI
379 <  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
380 < #else
381 <  ndfRaw = ndfRaw_local;
382 < #endif
383 <
384 <  return ndfRaw;
385 < }
386 <
387 < int SimInfo::getNDFtranslational() {
388 <  int ndfTrans_local;
389 <
390 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
391 <
392 <
393 < #ifdef IS_MPI
394 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
395 < #else
396 <  ndfTrans = ndfTrans_local;
397 < #endif
398 <
399 <  ndfTrans = ndfTrans - 3 - nZconstraints;
400 <
401 <  return ndfTrans;
402 < }
403 <
404 < int SimInfo::getTotIntegrableObjects() {
405 <  int nObjs_local;
406 <  int nObjs;
407 <
408 <  nObjs_local =  integrableObjects.size();
409 <
410 <
411 < #ifdef IS_MPI
412 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
413 < #else
414 <  nObjs = nObjs_local;
415 < #endif
416 <
417 <
418 <  return nObjs;
419 < }
420 <
421 < void SimInfo::refreshSim(){
422 <
423 <  simtype fInfo;
424 <  int isError;
425 <  int n_global;
426 <  int* excl;
427 <
428 <  fInfo.dielect = 0.0;
429 <
430 <  if( useDipoles ){
431 <    if( useReactionField )fInfo.dielect = dielectric;
432 <  }
433 <
434 <  fInfo.SIM_uses_PBC = usePBC;
435 <  //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;
446 <
447 <  n_exclude = excludes->getSize();
448 <  excl = excludes->getFortranArray();
449 <  
450 < #ifdef IS_MPI
451 <  n_global = mpiSim->getNAtomsGlobal();
452 < #else
453 <  n_global = n_atoms;
454 < #endif
455 <  
456 <  isError = 0;
457 <  
458 <  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
459 <  //it may not be a good idea to pass the address of first element in vector
460 <  //since c++ standard does not require vector to be stored continuously in meomory
461 <  //Most of the compilers will organize the memory of vector continuously
462 <  setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
463 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
464 <                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
465 <
466 <  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();
484 < }
485 <
486 < void SimInfo::setDefaultRcut( double theRcut ){
487 <  
488 <  haveRcut = 1;
489 <  rCut = theRcut;
490 <  rList = rCut + 1.0;
491 <  
492 <  notifyFortranCutOffs( &rCut, &rSw, &rList );
493 < }
494 <
495 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
496 <
497 <  rSw = theRsw;
498 <  setDefaultRcut( theRcut );
499 < }
500 <
501 <
502 < void SimInfo::checkCutOffs( void ){
503 <  
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 <  
536 < }
537 <
538 < void SimInfo::addProperty(GenericData* prop){
539 <
540 <  map<string, GenericData*>::iterator result;
541 <  result = properties.find(prop->getID());
542 <  
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;
550 <      
551 <  }
552 <  else{
553 <
554 <    properties[prop->getID()] = prop;
555 <
556 <  }
557 <    
558 < }
559 <
560 < GenericData* SimInfo::getProperty(const string& propName){
196 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
197  
562 –  map<string, GenericData*>::iterator result;
563 –  
564 –  //string lowerCaseName = ();
565 –  
566 –  result = properties.find(propName);
567 –  
568 –  if(result != properties.end())
569 –    return (*result).second;  
570 –  else  
571 –    return NULL;  
198   }
199  
200  
201 < void SimInfo::getFortranGroupArrays(SimInfo* info,
576 <                                    vector<int>& FglobalGroupMembership,
577 <                                    vector<double>& mfact){
578 <  
579 <  Molecule* myMols;
580 <  Atom** myAtoms;
581 <  int numAtom;
582 <  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 <  
595 <
596 <  // Fix the silly fortran indexing problem
597 < #ifdef IS_MPI
598 <  numAtom = mpiSim->getNAtomsGlobal();
599 < #else
600 <  numAtom = n_atoms;
601 < #endif
602 <  for (int i = 0; i < numAtom; i++)
603 <    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
604 <  
605 <
606 <  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)){
613 <
614 <      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 <  }
623 <
624 < }
201 > }//end namespace oopse

Diff Legend

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