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 branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents):
Revision 1683, Thu Oct 28 22:34:02 2004 UTC vs.
Revision 1719 by tim, Fri Nov 5 23:38:27 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  
33 – SimInfo::SimInfo(){
34 –
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;
48 –
49 –  haveRcut = 0;
50 –  haveRsw = 0;
51 –  boxIsInit = 0;
52 –  
53 –  resetTime = 1e99;
54 –
55 –  orthoRhombic = 0;
56 –  orthoTolerance = 1E-6;
57 –  useInitXSstate = true;
58 –
59 –  usePBC = 0;
60 –  useDirectionalAtoms = 0;
61 –  useLennardJones = 0;
62 –  useElectrostatics = 0;
63 –  useCharges = 0;
64 –  useDipoles = 0;
65 –  useSticky = 0;
66 –  useGayBerne = 0;
67 –  useEAM = 0;
68 –  useShapes = 0;
69 –  useFLARB = 0;
70 –
71 –  useSolidThermInt = 0;
72 –  useLiquidThermInt = 0;
73 –
74 –  haveCutoffGroups = false;
75 –
76 –  excludes = Exclude::Instance();
77 –
78 –  myConfiguration = new SimState();
79 –
80 –  has_minimizer = false;
81 –  the_minimizer =NULL;
82 –
83 –  ngroup = 0;
84 –
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;
96 <
67 >        return true;
68 >    } else {
69 >        return false;
70 >    }
71   }
72  
73 < void SimInfo::setBox(double newBox[3]) {
74 <  
75 <  int i, j;
102 <  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];
110 <
111 <  setBoxM( tempMat );
112 <
113 < }
114 <
115 < void SimInfo::setBoxM( double theBox[3][3] ){
116 <  
117 <  int i, j;
118 <  double FortranHmat[9]; // to preserve compatibility with Fortran the
119 <                         // ordering in the array is as follows:
120 <                         // [ 0 3 6 ]
121 <                         // [ 1 4 7 ]
122 <                         // [ 2 5 8 ]
123 <  double FortranHmatInv[9]; // the inverted Hmat (for Fortran);
124 <
125 <  if( !boxIsInit ) boxIsInit = 1;
126 <
127 <  for(i=0; i < 3; i++)
128 <    for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j];
129 <  
130 <  calcBoxL();
131 <  calcHmatInv();
132 <
133 <  for(i=0; i < 3; i++) {
134 <    for (j=0; j < 3; j++) {
135 <      FortranHmat[3*j + i] = Hmat[i][j];
136 <      FortranHmatInv[3*j + i] = HmatInv[i][j];
88 >        return true;
89 >    } else {
90 >        return false;
91      }
138 –  }
92  
140 –  setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic);
141 –
142 – }
143 –
93  
94 < void SimInfo::getBoxM (double theBox[3][3]) {
94 > }    
95  
96 <  int i, j;
97 <  for(i=0; i<3; i++)
98 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j];
99 < }
96 >        
97 > Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
98 >    i = molecules_.begin();
99 >    return i == molecules_.end() ? NULL : *i;
100 > }    
101  
102 <
103 < void SimInfo::scaleBox(double scale) {
104 <  double theBox[3][3];
155 <  int i, j;
156 <
157 <  // cerr << "Scaling box by " << scale << "\n";
158 <
159 <  for(i=0; i<3; i++)
160 <    for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale;
161 <
162 <  setBoxM(theBox);
163 <
102 > Molecule* SimInfo::nextMolecule(std::vector<Molecule*>::iterator& i) {
103 >    ++i;
104 >    return i == molecules_.end() ? NULL : *i;    
105   }
106  
166 – void SimInfo::calcHmatInv( void ) {
167 –  
168 –  int oldOrtho;
169 –  int i,j;
170 –  double smallDiag;
171 –  double tol;
172 –  double sanity[3][3];
107  
108 <  invertMat3( Hmat, HmatInv );
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 <  // check to see if Hmat is orthorhombic
177 <  
178 <  oldOrtho = orthoRhombic;
179 <
180 <  smallDiag = fabs(Hmat[0][0]);
181 <  if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]);
182 <  if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]);
183 <  tol = smallDiag * orthoTolerance;
184 <
185 <  orthoRhombic = 1;
186 <  
187 <  for (i = 0; i < 3; i++ ) {
188 <    for (j = 0 ; j < 3; j++) {
189 <      if (i != j) {
190 <        if (orthoRhombic) {
191 <          if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0;
192 <        }        
193 <      }
194 <    }
195 <  }
196 <
197 <  if( oldOrtho != orthoRhombic ){
115 >    ndf_local = 0;
116      
117 <    if( orthoRhombic ) {
118 <      sprintf( painCave.errMsg,
119 <               "OOPSE is switching from the default Non-Orthorhombic\n"
202 <               "\tto the faster Orthorhombic periodic boundary computations.\n"
203 <               "\tThis is usually a good thing, but if you wan't the\n"
204 <               "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n"
205 <               "\tvariable ( currently set to %G ) smaller.\n",
206 <               orthoTolerance);
207 <      painCave.severity = OOPSE_INFO;
208 <      simError();
209 <    }
210 <    else {
211 <      sprintf( painCave.errMsg,
212 <               "OOPSE is switching from the faster Orthorhombic to the more\n"
213 <               "\tflexible Non-Orthorhombic periodic boundary computations.\n"
214 <               "\tThis is usually because the box has deformed under\n"
215 <               "\tNPTf integration. If you wan't to live on the edge with\n"
216 <               "\tthe Orthorhombic computations, make the orthoBoxTolerance\n"
217 <               "\tvariable ( currently set to %G ) larger.\n",
218 <               orthoTolerance);
219 <      painCave.severity = OOPSE_WARNING;
220 <      simError();
221 <    }
222 <  }
223 < }
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::calcBoxL( void ){
121 >            ndf_local += 3;
122  
123 <  double dx, dy, dz, dsq;
124 <
125 <  // boxVol = Determinant of Hmat
126 <
127 <  boxVol = matDet3( Hmat );
128 <
129 <  // boxLx
130 <  
131 <  dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0];
132 <  dsq = dx*dx + dy*dy + dz*dz;
237 <  boxL[0] = sqrt( dsq );
238 <  //maxCutoff = 0.5 * boxL[0];
239 <
240 <  // boxLy
241 <  
242 <  dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1];
243 <  dsq = dx*dx + dy*dy + dz*dz;
244 <  boxL[1] = sqrt( dsq );
245 <  //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1];
246 <
247 <
248 <  // boxLz
249 <  
250 <  dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2];
251 <  dsq = dx*dx + dy*dy + dz*dz;
252 <  boxL[2] = sqrt( dsq );
253 <  //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2];
254 <
255 <  //calculate the max cutoff
256 <  maxCutoff =  calcMaxCutOff();
257 <  
258 <  checkCutOffs();
259 <
260 < }
261 <
262 <
263 < double SimInfo::calcMaxCutOff(){
264 <
265 <  double ri[3], rj[3], rk[3];
266 <  double rij[3], rjk[3], rki[3];
267 <  double minDist;
268 <
269 <  ri[0] = Hmat[0][0];
270 <  ri[1] = Hmat[1][0];
271 <  ri[2] = Hmat[2][0];
272 <
273 <  rj[0] = Hmat[0][1];
274 <  rj[1] = Hmat[1][1];
275 <  rj[2] = Hmat[2][1];
276 <
277 <  rk[0] = Hmat[0][2];
278 <  rk[1] = Hmat[1][2];
279 <  rk[2] = Hmat[2][2];
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 <  crossProduct3(ri, rj, rij);
135 <  distXY = dotProduct3(rk,rij) / norm3(rij);
134 >    // n_constraints is local, so subtract them on each processor
135 >    ndf_local -= nConstraints_;
136  
137 <  crossProduct3(rj,rk, rjk);
138 <  distYZ = dotProduct3(ri,rjk) / norm3(rjk);
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 <  crossProduct3(rk,ri, rki);
144 <  distZX = dotProduct3(rj,rki) / norm3(rki);
143 >    // nZconstraints is global, as are the 3 COM translations for the
144 >    // entire system:
145 >    ndf_ = ndf_ - 3 - nZconstraints;
146  
290 –  minDist = min(min(distXY, distYZ), distZX);
291 –  return minDist/2;
292 –  
147   }
148  
149 < void SimInfo::wrapVector( double thePos[3] ){
149 > void SimInfo::calcNDFRaw() {
150 >    int ndfRaw_local;
151  
152 <  int i;
153 <  double scaled[3];
152 >    std::vector<Molecule*>::iterator i;
153 >    std::vector<StuntDouble*>::iterator j;
154 >    Molecule* mol;
155 >    StuntDouble* integrableObject;
156  
157 <  if( !orthoRhombic ){
158 <    // calc the scaled coordinates.
302 <  
303 <
304 <    matVecMul3(HmatInv, thePos, scaled);
157 >    // Raw degrees of freedom that we have to set
158 >    ndfRaw_local = 0;
159      
160 <    for(i=0; i<3; i++)
161 <      scaled[i] -= roundMe(scaled[i]);
162 <    
309 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
310 <    
311 <    matVecMul3(Hmat, scaled, thePos);
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 <  }
314 <  else{
315 <    // calc the scaled coordinates.
316 <    
317 <    for(i=0; i<3; i++)
318 <      scaled[i] = thePos[i]*HmatInv[i][i];
319 <    
320 <    // wrap the scaled coordinates
321 <    
322 <    for(i=0; i<3; i++)
323 <      scaled[i] -= roundMe(scaled[i]);
324 <    
325 <    // calc the wrapped real coordinates from the wrapped scaled coordinates
326 <    
327 <    for(i=0; i<3; i++)
328 <      thePos[i] = scaled[i]*Hmat[i][i];
329 <  }
330 <    
331 < }
164 >            ndfRaw_local += 3;
165  
166 <
167 < int SimInfo::getNDF(){
168 <  int ndf_local;
169 <
170 <  ndf_local = 0;
171 <  
172 <  for(int i = 0; i < integrableObjects.size(); i++){
173 <    ndf_local += 3;
174 <    if (integrableObjects[i]->isDirectional()) {
342 <      if (integrableObjects[i]->isLinear())
343 <        ndf_local += 2;
344 <      else
345 <        ndf_local += 3;
166 >            if (integrableObject->isDirectional()) {
167 >                if (integrableObject->isLinear()) {
168 >                    ndfRaw_local += 2;
169 >                } else {
170 >                    ndfRaw_local += 3;
171 >                }
172 >            }
173 >            
174 >        }
175      }
347 –  }
348 –
349 –  // n_constraints is local, so subtract them on each processor:
350 –
351 –  ndf_local -= n_constraints;
352 –
353 – #ifdef IS_MPI
354 –  MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
355 – #else
356 –  ndf = ndf_local;
357 – #endif
358 –
359 –  // nZconstraints is global, as are the 3 COM translations for the
360 –  // entire system:
361 –
362 –  ndf = ndf - 3 - nZconstraints;
363 –
364 –  return ndf;
365 – }
366 –
367 – int SimInfo::getNDFraw() {
368 –  int ndfRaw_local;
369 –
370 –  // Raw degrees of freedom that we have to set
371 –  ndfRaw_local = 0;
372 –
373 –  for(int i = 0; i < integrableObjects.size(); i++){
374 –    ndfRaw_local += 3;
375 –    if (integrableObjects[i]->isDirectional()) {
376 –       if (integrableObjects[i]->isLinear())
377 –        ndfRaw_local += 2;
378 –      else
379 –        ndfRaw_local += 3;
380 –    }
381 –  }
176      
177   #ifdef IS_MPI
178 <  MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
178 >    MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
179   #else
180 <  ndfRaw = ndfRaw_local;
180 >    ndfRaw_ = ndfRaw_local;
181   #endif
388 –
389 –  return ndfRaw;
182   }
183  
184 < int SimInfo::getNDFtranslational() {
185 <  int ndfTrans_local;
184 > void SimInfo::calcNDFTrans() {
185 >    int ndfTrans_local;
186  
187 <  ndfTrans_local = 3 * integrableObjects.size() - n_constraints;
187 >    ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
188  
189  
190   #ifdef IS_MPI
191 <  MPI_Allreduce(&ndfTrans_local,&ndfTrans,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 <  ndfTrans = ndfTrans_local;
193 >    ndfTrans_ = ndfTrans_local;
194   #endif
195  
196 <  ndfTrans = ndfTrans - 3 - nZconstraints;
197 <
406 <  return ndfTrans;
196 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
197 >
198   }
199  
200 < int SimInfo::getTotIntegrableObjects() {
201 <  int nObjs_local;
202 <  int nObjs;
203 <
204 <  nObjs_local =  integrableObjects.size();
205 <
206 <
207 < #ifdef IS_MPI
208 <  MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
209 < #else
210 <  nObjs = nObjs_local;
420 < #endif
421 <
422 <
423 <  return nObjs;
424 < }
425 <
426 < void SimInfo::refreshSim(){
427 <
428 <  simtype fInfo;
429 <  int isError;
430 <  int n_global;
431 <  int* excl;
432 <
433 <  fInfo.dielect = 0.0;
434 <
435 <  if( useDipoles ){
436 <    if( useReactionField )fInfo.dielect = dielectric;
437 <  }
438 <
439 <  fInfo.SIM_uses_PBC = usePBC;
440 <
441 <  if (useSticky || useDipoles || useGayBerne || useShapes) {
442 <    useDirectionalAtoms = 1;
443 <    fInfo.SIM_uses_DirectionalAtoms = useDirectionalAtoms;
444 <  }
445 <
446 <  fInfo.SIM_uses_LennardJones = useLennardJones;
447 <
448 <  if (useCharges || useDipoles) {
449 <    useElectrostatics = 1;
450 <    fInfo.SIM_uses_Electrostatics = useElectrostatics;
451 <  }
452 <
453 <  fInfo.SIM_uses_Charges = useCharges;
454 <  fInfo.SIM_uses_Dipoles = useDipoles;
455 <  fInfo.SIM_uses_Sticky = useSticky;
456 <  fInfo.SIM_uses_GayBerne = useGayBerne;
457 <  fInfo.SIM_uses_EAM = useEAM;
458 <  fInfo.SIM_uses_Shapes = useShapes;
459 <  fInfo.SIM_uses_FLARB = useFLARB;
460 <  fInfo.SIM_uses_RF = useReactionField;
461 <
462 <  n_exclude = excludes->getSize();
463 <  excl = excludes->getFortranArray();
464 <  
465 < #ifdef IS_MPI
466 <  n_global = mpiSim->getNAtomsGlobal();
467 < #else
468 <  n_global = n_atoms;
469 < #endif
470 <  
471 <  isError = 0;
472 <  
473 <  getFortranGroupArrays(this, FglobalGroupMembership, mfact);
474 <  //it may not be a good idea to pass the address of first element in vector
475 <  //since c++ standard does not require vector to be stored continuously in meomory
476 <  //Most of the compilers will organize the memory of vector continuously
477 <  setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl,
478 <                  &nGlobalExcludes, globalExcludes, molMembershipArray,
479 <                  &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError);
480 <
481 <  if( isError ){
200 > void SimInfo::addExcludePairs(Molecule* mol) {
201 >    std::vector<Bond*>::iterator bondIter;
202 >    std::vector<Bend*>::iterator bendIter;
203 >    std::vector<Torsion*>::iterator torsionIter;
204 >    Bond* bond;
205 >    Bend* bend;
206 >    Torsion* torsion;
207 >    int a;
208 >    int b;
209 >    int c;
210 >    int d;
211      
212 <    sprintf( painCave.errMsg,
213 <             "There was an error setting the simulation information in fortran.\n" );
214 <    painCave.isFatal = 1;
215 <    painCave.severity = OOPSE_ERROR;
216 <    simError();
488 <  }
489 <  
490 < #ifdef IS_MPI
491 <  sprintf( checkPointMsg,
492 <           "succesfully sent the simulation information to fortran.\n");
493 <  MPIcheckPoint();
494 < #endif // is_mpi
495 <  
496 <  this->ndf = this->getNDF();
497 <  this->ndfRaw = this->getNDFraw();
498 <  this->ndfTrans = this->getNDFtranslational();
499 < }
212 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
213 >        a = bond->getAtomA()->getGlobalIndex();
214 >        b = bond->getAtomB()->getGlobalIndex();        
215 >        exclude_.addPair(a, b);
216 >    }
217  
218 < void SimInfo::setDefaultRcut( double theRcut ){
219 <  
220 <  haveRcut = 1;
221 <  rCut = theRcut;
505 <  rList = rCut + 1.0;
506 <  
507 <  notifyFortranCutoffs( &rCut, &rSw, &rList );
508 < }
218 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
219 >        a = bend->getAtomA()->getGlobalIndex();
220 >        b = bend->getAtomB()->getGlobalIndex();        
221 >        c = bend->getAtomC()->getGlobalIndex();
222  
223 < void SimInfo::setDefaultRcut( double theRcut, double theRsw ){
223 >        exclude_.addPair(a, b);
224 >        exclude_.addPair(a, c);
225 >        exclude_.addPair(b, c);        
226 >    }
227  
228 <  rSw = theRsw;
229 <  setDefaultRcut( theRcut );
230 < }
228 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
229 >        a = torsion->getAtomA()->getGlobalIndex();
230 >        b = torsion->getAtomB()->getGlobalIndex();        
231 >        c = torsion->getAtomC()->getGlobalIndex();        
232 >        d = torsion->getAtomD()->getGlobalIndex();        
233  
234 +        exclude_.addPair(a, b);
235 +        exclude_.addPair(a, c);
236 +        exclude_.addPair(a, d);
237 +        exclude_.addPair(b, c);
238 +        exclude_.addPair(b, d);
239 +        exclude_.addPair(c, d);        
240 +    }
241  
517 – void SimInfo::checkCutOffs( void ){
518 –  
519 –  if( boxIsInit ){
242      
521 –    //we need to check cutOffs against the box
522 –    
523 –    if( rCut > maxCutoff ){
524 –      sprintf( painCave.errMsg,
525 –               "cutoffRadius is too large for the current periodic box.\n"
526 –               "\tCurrent Value of cutoffRadius = %G at time %G\n "
527 –               "\tThis is larger than half of at least one of the\n"
528 –               "\tperiodic box vectors.  Right now, the Box matrix is:\n"
529 –               "\n"
530 –               "\t[ %G %G %G ]\n"
531 –               "\t[ %G %G %G ]\n"
532 –               "\t[ %G %G %G ]\n",
533 –               rCut, currentTime,
534 –               Hmat[0][0], Hmat[0][1], Hmat[0][2],
535 –               Hmat[1][0], Hmat[1][1], Hmat[1][2],
536 –               Hmat[2][0], Hmat[2][1], Hmat[2][2]);
537 –      painCave.severity = OOPSE_ERROR;
538 –      painCave.isFatal = 1;
539 –      simError();
540 –    }    
541 –  } else {
542 –    // initialize this stuff before using it, OK?
543 –    sprintf( painCave.errMsg,
544 –             "Trying to check cutoffs without a box.\n"
545 –             "\tOOPSE should have better programmers than that.\n" );
546 –    painCave.severity = OOPSE_ERROR;
547 –    painCave.isFatal = 1;
548 –    simError();      
549 –  }
550 –  
243   }
244  
245 < void SimInfo::addProperty(GenericData* prop){
246 <
247 <  map<string, GenericData*>::iterator result;
248 <  result = properties.find(prop->getID());
249 <  
250 <  //we can't simply use  properties[prop->getID()] = prop,
251 <  //it will cause memory leak if we already contain a propery which has the same name of prop
252 <  
253 <  if(result != properties.end()){
245 > void SimInfo::removeExcludePairs(Molecule* mol) {
246 >    std::vector<Bond*>::iterator bondIter;
247 >    std::vector<Bend*>::iterator bendIter;
248 >    std::vector<Torsion*>::iterator torsionIter;
249 >    Bond* bond;
250 >    Bend* bend;
251 >    Torsion* torsion;
252 >    int a;
253 >    int b;
254 >    int c;
255 >    int d;
256      
257 <    delete (*result).second;
258 <    (*result).second = prop;
259 <      
260 <  }
261 <  else{
257 >    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
258 >        a = bond->getAtomA()->getGlobalIndex();
259 >        b = bond->getAtomB()->getGlobalIndex();        
260 >        exclude_.removePair(a, b);
261 >    }
262  
263 <    properties[prop->getID()] = prop;
263 >    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
264 >        a = bend->getAtomA()->getGlobalIndex();
265 >        b = bend->getAtomB()->getGlobalIndex();        
266 >        c = bend->getAtomC()->getGlobalIndex();
267  
268 <  }
269 <    
270 < }
268 >        exclude_.removePair(a, b);
269 >        exclude_.removePair(a, c);
270 >        exclude_.removePair(b, c);        
271 >    }
272  
273 < GenericData* SimInfo::getProperty(const string& propName){
274 <
275 <  map<string, GenericData*>::iterator result;
276 <  
277 <  //string lowerCaseName = ();
580 <  
581 <  result = properties.find(propName);
582 <  
583 <  if(result != properties.end())
584 <    return (*result).second;  
585 <  else  
586 <    return NULL;  
587 < }
273 >    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextBond(torsionIter)) {
274 >        a = torsion->getAtomA()->getGlobalIndex();
275 >        b = torsion->getAtomB()->getGlobalIndex();        
276 >        c = torsion->getAtomC()->getGlobalIndex();        
277 >        d = torsion->getAtomD()->getGlobalIndex();        
278  
279 <
280 < void SimInfo::getFortranGroupArrays(SimInfo* info,
281 <                                    vector<int>& FglobalGroupMembership,
282 <                                    vector<double>& mfact){
283 <  
284 <  Molecule* myMols;
595 <  Atom** myAtoms;
596 <  int numAtom;
597 <  double mtot;
598 <  int numMol;
599 <  int numCutoffGroups;
600 <  CutoffGroup* myCutoffGroup;
601 <  vector<CutoffGroup*>::iterator iterCutoff;
602 <  Atom* cutoffAtom;
603 <  vector<Atom*>::iterator iterAtom;
604 <  int atomIndex;
605 <  double totalMass;
606 <  
607 <  mfact.clear();
608 <  FglobalGroupMembership.clear();
609 <  
610 <
611 <  // Fix the silly fortran indexing problem
612 < #ifdef IS_MPI
613 <  numAtom = mpiSim->getNAtomsGlobal();
614 < #else
615 <  numAtom = n_atoms;
616 < #endif
617 <  for (int i = 0; i < numAtom; i++)
618 <    FglobalGroupMembership.push_back(globalGroupMembership[i] + 1);
619 <  
620 <
621 <  myMols = info->molecules;
622 <  numMol = info->n_mol;
623 <  for(int i  = 0; i < numMol; i++){
624 <    numCutoffGroups = myMols[i].getNCutoffGroups();
625 <    for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff);
626 <        myCutoffGroup != NULL;
627 <        myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){
628 <
629 <      totalMass = myCutoffGroup->getMass();
630 <      
631 <      for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom);
632 <          cutoffAtom != NULL;
633 <          cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){
634 <        mfact.push_back(cutoffAtom->getMass()/totalMass);
635 <      }  
279 >        exclude_.removePair(a, b);
280 >        exclude_.removePair(a, c);
281 >        exclude_.removePair(a, d);
282 >        exclude_.removePair(b, c);
283 >        exclude_.removePair(b, d);
284 >        exclude_.removePair(c, d);        
285      }
637 –  }
286  
287   }
288 +
289 +
290 + }//end namespace oopse

Diff Legend

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