ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/SimSetup.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/SimSetup.cpp (file contents):
Revision 708 by tim, Wed Aug 20 22:23:34 2003 UTC vs.
Revision 722 by tim, Tue Aug 26 20:02:29 2003 UTC

# Line 3 | Line 3
3   #include <iostream>
4   #include <cmath>
5   #include <string>
6 < #include <sys/time.h>
6 > #include <sprng.h>
7  
8   #include "SimSetup.hpp"
9   #include "ReadWrite.hpp"
# Line 32 | Line 32 | SimSetup::SimSetup(){
32   using namespace std;
33  
34   SimSetup::SimSetup(){
35  
35    isInfoArray = 0;
36    nInfo = 1;
37 <  
37 >
38    stamps = new MakeStamps();
39    globals = new Globals();
40 <  
41 <  
40 >
41 >
42   #ifdef IS_MPI
43 <  strcpy( checkPointMsg, "SimSetup creation successful" );
43 >  strcpy(checkPointMsg, "SimSetup creation successful");
44    MPIcheckPoint();
45   #endif // IS_MPI
46   }
# Line 51 | Line 50 | void SimSetup::setSimInfo( SimInfo* the_info, int theN
50    delete globals;
51   }
52  
53 < void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
54 <    info = the_info;
55 <    nInfo = theNinfo;
56 <    isInfoArray = 1;
53 > void SimSetup::setSimInfo(SimInfo* the_info, int theNinfo){
54 >  info = the_info;
55 >  nInfo = theNinfo;
56 >  isInfoArray = 1;
57   }
58  
59  
60 < void SimSetup::parseFile( char* fileName ){
62 <
60 > void SimSetup::parseFile(char* fileName){
61   #ifdef IS_MPI
62 <  if( worldRank == 0 ){
62 >  if (worldRank == 0){
63   #endif // is_mpi
64 <    
64 >
65      inFileName = fileName;
66 <    set_interface_stamps( stamps, globals );
67 <    
66 >    set_interface_stamps(stamps, globals);
67 >
68   #ifdef IS_MPI
69      mpiEventInit();
70   #endif
71  
72 <    yacc_BASS( fileName );
72 >    yacc_BASS(fileName);
73  
74   #ifdef IS_MPI
75      throwMPIEvent(NULL);
76    }
77 <  else receiveParse();
77 >  else{
78 >    receiveParse();
79 >  }
80   #endif
81  
82   }
83  
84   #ifdef IS_MPI
85   void SimSetup::receiveParse(void){
86 <
87 <    set_interface_stamps( stamps, globals );
88 <    mpiEventInit();
89 <    MPIcheckPoint();
90 <    mpiEventLoop();
91 <
86 >  set_interface_stamps(stamps, globals);
87 >  mpiEventInit();
88 >  MPIcheckPoint();
89 >  mpiEventLoop();
90   }
91  
92   #endif // is_mpi
93  
94   void SimSetup::createSim(void){
97
95    int i, j, k, globalAtomIndex;
96 <  
96 >
97    // gather all of the information from the Bass file
98  
99    gatherInfo();
# Line 111 | Line 108 | void SimSetup::createSim(void){
108  
109    // initialize the system coordinates
110  
111 <  if( !isInfoArray ) initSystemCoords();  
111 >  if (!isInfoArray){
112 >    initSystemCoords();
113 >  }  
114  
115    // make the output filenames
116  
117    makeOutNames();
118 <  
118 >
119    // make the integrator
120 <  
120 >
121    makeIntegrator();
122 <  
122 >
123   #ifdef IS_MPI
124    mpiSim->mpiRefresh();
125   #endif
# Line 128 | Line 127 | void SimSetup::createSim(void){
127    // initialize the Fortran
128  
129    initFortran();
131
132
133
130   }
131  
132  
133 < void SimSetup::makeMolecules( void ){
134 <
139 <  int k,l;
133 > void SimSetup::makeMolecules(void){
134 >  int k, l;
135    int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
136    molInit molInfo;
137    DirectionalAtom* dAtom;
# Line 151 | Line 146 | void SimSetup::makeMolecules( void ){
146    bend_set* theBends;
147    torsion_set* theTorsions;
148  
149 <  
149 >
150    //init the forceField paramters
151  
152    the_ff->readParams();
153  
154 <  
154 >
155    // init the atoms
156  
157    double ux, uy, uz, u, uSqr;
163  
164  for(k=0; k<nInfo; k++){
165    
166    the_ff->setSimInfo( &(info[k]) );
158  
159 +  for (k = 0; k < nInfo; k++){
160 +    the_ff->setSimInfo(&(info[k]));
161 +
162      atomOffset = 0;
163      excludeOffset = 0;
164 <    for(i=0; i<info[k].n_mol; i++){
171 <    
164 >    for (i = 0; i < info[k].n_mol; i++){
165        stampID = info[k].molecules[i].getStampID();
166  
167 <      molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
168 <      molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
169 <      molInfo.nBends    = comp_stamps[stampID]->getNBends();
167 >      molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
168 >      molInfo.nBonds = comp_stamps[stampID]->getNBonds();
169 >      molInfo.nBends = comp_stamps[stampID]->getNBends();
170        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
171        molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
172 <      
172 >
173        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
174        molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
175 <      molInfo.myBonds = new Bond*[molInfo.nBonds];
176 <      molInfo.myBends = new Bend*[molInfo.nBends];
177 <      molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
175 >      molInfo.myBonds = new Bond * [molInfo.nBonds];
176 >      molInfo.myBends = new Bend * [molInfo.nBends];
177 >      molInfo.myTorsions = new Torsion * [molInfo.nTorsions];
178  
179        theBonds = new bond_pair[molInfo.nBonds];
180        theBends = new bend_set[molInfo.nBends];
181        theTorsions = new torsion_set[molInfo.nTorsions];
182 <    
182 >
183        // make the Atoms
184 <    
185 <      for(j=0; j<molInfo.nAtoms; j++){
186 <  
187 <  currentAtom = comp_stamps[stampID]->getAtom( j );
188 <  if( currentAtom->haveOrientation() ){
189 <    
190 <    dAtom = new DirectionalAtom( (j + atomOffset),
191 <               info[k].getConfiguration() );
192 <    info[k].n_oriented++;
193 <    molInfo.myAtoms[j] = dAtom;
194 <    
195 <    ux = currentAtom->getOrntX();
196 <    uy = currentAtom->getOrntY();
197 <    uz = currentAtom->getOrntZ();
198 <    
199 <    uSqr = (ux * ux) + (uy * uy) + (uz * uz);
200 <    
201 <    u = sqrt( uSqr );
202 <    ux = ux / u;
203 <    uy = uy / u;
204 <    uz = uz / u;
205 <    
206 <    dAtom->setSUx( ux );
207 <    dAtom->setSUy( uy );
208 <    dAtom->setSUz( uz );
209 <  }
210 <  else{
211 <    molInfo.myAtoms[j] = new GeneralAtom( (j + atomOffset),
212 <            info[k].getConfiguration() );
213 <  }
221 <  molInfo.myAtoms[j]->setType( currentAtom->getType() );
222 <    
184 >
185 >      for (j = 0; j < molInfo.nAtoms; j++){
186 >        currentAtom = comp_stamps[stampID]->getAtom(j);
187 >        if (currentAtom->haveOrientation()){
188 >          dAtom = new DirectionalAtom((j + atomOffset),
189 >                                      info[k].getConfiguration());
190 >          info[k].n_oriented++;
191 >          molInfo.myAtoms[j] = dAtom;
192 >
193 >          ux = currentAtom->getOrntX();
194 >          uy = currentAtom->getOrntY();
195 >          uz = currentAtom->getOrntZ();
196 >
197 >          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
198 >
199 >          u = sqrt(uSqr);
200 >          ux = ux / u;
201 >          uy = uy / u;
202 >          uz = uz / u;
203 >
204 >          dAtom->setSUx(ux);
205 >          dAtom->setSUy(uy);
206 >          dAtom->setSUz(uz);
207 >        }
208 >        else{
209 >          molInfo.myAtoms[j] = new GeneralAtom((j + atomOffset),
210 >                                               info[k].getConfiguration());
211 >        }
212 >        molInfo.myAtoms[j]->setType(currentAtom->getType());
213 >
214   #ifdef IS_MPI
215 <      
216 <  molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
217 <      
215 >
216 >        molInfo.myAtoms[j]->setGlobalIndex(globalIndex[j + atomOffset]);
217 >
218   #endif // is_mpi
219        }
220 <    
221 <    // make the bonds
222 <      for(j=0; j<molInfo.nBonds; j++){
223 <      
224 <  currentBond = comp_stamps[stampID]->getBond( j );
225 <  theBonds[j].a = currentBond->getA() + atomOffset;
226 <  theBonds[j].b = currentBond->getB() + atomOffset;
227 <  
228 <  exI = theBonds[j].a;
229 <  exJ = theBonds[j].b;
230 <  
231 <  // exclude_I must always be the smaller of the pair
232 <  if( exI > exJ ){
233 <    tempEx = exI;
234 <    exI = exJ;
235 <    exJ = tempEx;
245 <  }
220 >
221 >      // make the bonds
222 >      for (j = 0; j < molInfo.nBonds; j++){
223 >        currentBond = comp_stamps[stampID]->getBond(j);
224 >        theBonds[j].a = currentBond->getA() + atomOffset;
225 >        theBonds[j].b = currentBond->getB() + atomOffset;
226 >
227 >        exI = theBonds[j].a;
228 >        exJ = theBonds[j].b;
229 >
230 >        // exclude_I must always be the smaller of the pair
231 >        if (exI > exJ){
232 >          tempEx = exI;
233 >          exI = exJ;
234 >          exJ = tempEx;
235 >        }
236   #ifdef IS_MPI
237 <  tempEx = exI;
238 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
239 <  tempEx = exJ;
240 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
241 <  
242 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
237 >        tempEx = exI;
238 >        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
239 >        tempEx = exJ;
240 >        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
241 >
242 >        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
243   #else  // isn't MPI
244 <  
245 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
244 >
245 >        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
246   #endif  //is_mpi
247        }
248        excludeOffset += molInfo.nBonds;
249 <      
249 >
250        //make the bends
251 <      for(j=0; j<molInfo.nBends; j++){
252 <  
253 <  currentBend = comp_stamps[stampID]->getBend( j );
254 <  theBends[j].a = currentBend->getA() + atomOffset;
255 <  theBends[j].b = currentBend->getB() + atomOffset;
256 <  theBends[j].c = currentBend->getC() + atomOffset;
257 <  
258 <  if( currentBend->haveExtras() ){
259 <    
260 <    extras = currentBend->getExtras();
261 <    current_extra = extras;
262 <    
263 <    while( current_extra != NULL ){
264 <      if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
265 <        
266 <        switch( current_extra->getType() ){
267 <    
268 <        case 0:
269 <    theBends[j].ghost =
270 <      current_extra->getInt() + atomOffset;
271 <    theBends[j].isGhost = 1;
272 <    break;
273 <    
274 <        case 1:
275 <    theBends[j].ghost =
276 <      (int)current_extra->getDouble() + atomOffset;
277 <    theBends[j].isGhost = 1;
278 <    break;
279 <    
280 <        default:
281 <    sprintf( painCave.errMsg,
282 <       "SimSetup Error: ghostVectorSource was neither a "
283 <       "double nor an int.\n"
284 <       "-->Bend[%d] in %s\n",
285 <       j, comp_stamps[stampID]->getID() );
286 <    painCave.isFatal = 1;
287 <    simError();
288 <        }
289 <      }
290 <      
291 <      else{
292 <        
293 <        sprintf( painCave.errMsg,
294 <           "SimSetup Error: unhandled bend assignment:\n"
295 <           "    -->%s in Bend[%d] in %s\n",
296 <           current_extra->getlhs(),
297 <           j, comp_stamps[stampID]->getID() );
298 <        painCave.isFatal = 1;
299 <        simError();
300 <      }
301 <      
302 <      current_extra = current_extra->getNext();
303 <    }
304 <  }
305 <  
306 <  if( !theBends[j].isGhost ){
307 <    
308 <    exI = theBends[j].a;
309 <    exJ = theBends[j].c;
310 <  }
311 <  else{
312 <    
323 <    exI = theBends[j].a;
324 <    exJ = theBends[j].b;
325 <  }
326 <  
327 <  // exclude_I must always be the smaller of the pair
328 <  if( exI > exJ ){
329 <    tempEx = exI;
330 <    exI = exJ;
331 <    exJ = tempEx;
332 <  }
251 >      for (j = 0; j < molInfo.nBends; j++){
252 >        currentBend = comp_stamps[stampID]->getBend(j);
253 >        theBends[j].a = currentBend->getA() + atomOffset;
254 >        theBends[j].b = currentBend->getB() + atomOffset;
255 >        theBends[j].c = currentBend->getC() + atomOffset;
256 >
257 >        if (currentBend->haveExtras()){
258 >          extras = currentBend->getExtras();
259 >          current_extra = extras;
260 >
261 >          while (current_extra != NULL){
262 >            if (!strcmp(current_extra->getlhs(), "ghostVectorSource")){
263 >              switch (current_extra->getType()){
264 >                case 0:
265 >                  theBends[j].ghost = current_extra->getInt() + atomOffset;
266 >                  theBends[j].isGhost = 1;
267 >                  break;
268 >
269 >                case 1:
270 >                  theBends[j].ghost = (int) current_extra->getDouble() +
271 >                                      atomOffset;
272 >                  theBends[j].isGhost = 1;
273 >                  break;
274 >
275 >                default:
276 >                  sprintf(painCave.errMsg,
277 >                          "SimSetup Error: ghostVectorSource was neither a "
278 >                          "double nor an int.\n"
279 >                          "-->Bend[%d] in %s\n",
280 >                          j, comp_stamps[stampID]->getID());
281 >                  painCave.isFatal = 1;
282 >                  simError();
283 >              }
284 >            }
285 >            else{
286 >              sprintf(painCave.errMsg,
287 >                      "SimSetup Error: unhandled bend assignment:\n"
288 >                      "    -->%s in Bend[%d] in %s\n",
289 >                      current_extra->getlhs(), j, comp_stamps[stampID]->getID());
290 >              painCave.isFatal = 1;
291 >              simError();
292 >            }
293 >
294 >            current_extra = current_extra->getNext();
295 >          }
296 >        }
297 >
298 >        if (!theBends[j].isGhost){
299 >          exI = theBends[j].a;
300 >          exJ = theBends[j].c;
301 >        }
302 >        else{
303 >          exI = theBends[j].a;
304 >          exJ = theBends[j].b;
305 >        }
306 >
307 >        // exclude_I must always be the smaller of the pair
308 >        if (exI > exJ){
309 >          tempEx = exI;
310 >          exI = exJ;
311 >          exJ = tempEx;
312 >        }
313   #ifdef IS_MPI
314 <  tempEx = exI;
315 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
316 <  tempEx = exJ;
317 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
318 <      
319 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
314 >        tempEx = exI;
315 >        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
316 >        tempEx = exJ;
317 >        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
318 >
319 >        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
320   #else  // isn't MPI
321 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
321 >        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
322   #endif  //is_mpi
323        }
324        excludeOffset += molInfo.nBends;
325 <      
326 <      for(j=0; j<molInfo.nTorsions; j++){
327 <  
328 <  currentTorsion = comp_stamps[stampID]->getTorsion( j );
329 <  theTorsions[j].a = currentTorsion->getA() + atomOffset;
330 <  theTorsions[j].b = currentTorsion->getB() + atomOffset;
331 <  theTorsions[j].c = currentTorsion->getC() + atomOffset;
332 <  theTorsions[j].d = currentTorsion->getD() + atomOffset;
333 <  
334 <  exI = theTorsions[j].a;
335 <  exJ = theTorsions[j].d;
336 <  
337 <  // exclude_I must always be the smaller of the pair
338 <  if( exI > exJ ){
339 <    tempEx = exI;
340 <    exI = exJ;
341 <    exJ = tempEx;
362 <  }
325 >
326 >      for (j = 0; j < molInfo.nTorsions; j++){
327 >        currentTorsion = comp_stamps[stampID]->getTorsion(j);
328 >        theTorsions[j].a = currentTorsion->getA() + atomOffset;
329 >        theTorsions[j].b = currentTorsion->getB() + atomOffset;
330 >        theTorsions[j].c = currentTorsion->getC() + atomOffset;
331 >        theTorsions[j].d = currentTorsion->getD() + atomOffset;
332 >
333 >        exI = theTorsions[j].a;
334 >        exJ = theTorsions[j].d;
335 >
336 >        // exclude_I must always be the smaller of the pair
337 >        if (exI > exJ){
338 >          tempEx = exI;
339 >          exI = exJ;
340 >          exJ = tempEx;
341 >        }
342   #ifdef IS_MPI
343 <  tempEx = exI;
344 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
345 <  tempEx = exJ;
346 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
347 <  
348 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
343 >        tempEx = exI;
344 >        exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
345 >        tempEx = exJ;
346 >        exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
347 >
348 >        info[k].excludes[j + excludeOffset]->setPair(exI, exJ);
349   #else  // isn't MPI
350 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
350 >        info[k].excludes[j + excludeOffset]->setPair((exI + 1), (exJ + 1));
351   #endif  //is_mpi
352        }
353        excludeOffset += molInfo.nTorsions;
354 <      
355 <      
354 >
355 >
356        // send the arrays off to the forceField for init.
378      
379      the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
380      the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
381      the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
382      the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
383      
384      
385      info[k].molecules[i].initialize( molInfo );
357  
358 <      
358 >      the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
359 >      the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
360 >      the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
361 >      the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
362 >                                 theTorsions);
363 >
364 >
365 >      info[k].molecules[i].initialize(molInfo);
366 >
367 >
368        atomOffset += molInfo.nAtoms;
369        delete[] theBonds;
370        delete[] theBends;
371        delete[] theTorsions;
372      }
373    }
374 <  
374 >
375   #ifdef IS_MPI
376 <  sprintf( checkPointMsg, "all molecules initialized succesfully" );
376 >  sprintf(checkPointMsg, "all molecules initialized succesfully");
377    MPIcheckPoint();
378   #endif // is_mpi
379 <  
379 >
380    // clean up the forcefield
381  
382    the_ff->calcRcut();
383    the_ff->cleanMe();
404  
384   }
385  
386 < void SimSetup::initFromBass( void ){
408 <
386 > void SimSetup::initFromBass(void){
387    int i, j, k;
388    int n_cells;
389    double cellx, celly, cellz;
# Line 419 | Line 397 | void SimSetup::initFromBass( void ){
397    vel[1] = 0.0;
398    vel[2] = 0.0;
399  
400 <  temp1 = (double)tot_nmol / 4.0;
401 <  temp2 = pow( temp1, ( 1.0 / 3.0 ) );
402 <  temp3 = ceil( temp2 );
400 >  temp1 = (double) tot_nmol / 4.0;
401 >  temp2 = pow(temp1, (1.0 / 3.0));
402 >  temp3 = ceil(temp2);
403  
404 <  have_extra =0;
405 <  if( temp2 < temp3 ){ // we have a non-complete lattice
406 <    have_extra =1;
404 >  have_extra = 0;
405 >  if (temp2 < temp3){
406 >    // we have a non-complete lattice
407 >    have_extra = 1;
408  
409 <    n_cells = (int)temp3 - 1;
409 >    n_cells = (int) temp3 - 1;
410      cellx = info[0].boxL[0] / temp3;
411      celly = info[0].boxL[1] / temp3;
412      cellz = info[0].boxL[2] / temp3;
413 <    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
414 <    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
415 <    n_per_extra = (int)ceil( temp1 );
413 >    n_extra = tot_nmol - (4 * n_cells * n_cells * n_cells);
414 >    temp1 = ((double) n_extra) / (pow(temp3, 3.0) - pow(n_cells, 3.0));
415 >    n_per_extra = (int) ceil(temp1);
416  
417 <    if( n_per_extra > 4){
418 <      sprintf( painCave.errMsg,
419 <         "SimSetup error. There has been an error in constructing"
420 <         " the non-complete lattice.\n" );
417 >    if (n_per_extra > 4){
418 >      sprintf(painCave.errMsg,
419 >              "SimSetup error. There has been an error in constructing"
420 >              " the non-complete lattice.\n");
421        painCave.isFatal = 1;
422        simError();
423      }
424    }
425    else{
426 <    n_cells = (int)temp3;
426 >    n_cells = (int) temp3;
427      cellx = info[0].boxL[0] / temp3;
428      celly = info[0].boxL[1] / temp3;
429      cellz = info[0].boxL[2] / temp3;
# Line 455 | Line 434 | void SimSetup::initFromBass( void ){
434    current_comp = 0;
435    current_atom_ndx = 0;
436  
437 <  for( i=0; i < n_cells ; i++ ){
438 <    for( j=0; j < n_cells; j++ ){
439 <      for( k=0; k < n_cells; k++ ){
437 >  for (i = 0; i < n_cells ; i++){
438 >    for (j = 0; j < n_cells; j++){
439 >      for (k = 0; k < n_cells; k++){
440 >        makeElement(i * cellx, j * celly, k * cellz);
441  
442 <  makeElement( i * cellx,
463 <         j * celly,
464 <         k * cellz );
442 >        makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly, k * cellz);
443  
444 <  makeElement( i * cellx + 0.5 * cellx,
467 <         j * celly + 0.5 * celly,
468 <         k * cellz );
444 >        makeElement(i * cellx, j * celly + 0.5 * celly, k * cellz + 0.5 * cellz);
445  
446 <  makeElement( i * cellx,
471 <         j * celly + 0.5 * celly,
472 <         k * cellz + 0.5 * cellz );
473 <
474 <  makeElement( i * cellx + 0.5 * cellx,
475 <         j * celly,
476 <         k * cellz + 0.5 * cellz );
446 >        makeElement(i * cellx + 0.5 * cellx, j * celly, k * cellz + 0.5 * cellz);
447        }
448      }
449    }
450  
451 <  if( have_extra ){
451 >  if (have_extra){
452      done = 0;
453  
454      int start_ndx;
455 <    for( i=0; i < (n_cells+1) && !done; i++ ){
456 <      for( j=0; j < (n_cells+1) && !done; j++ ){
455 >    for (i = 0; i < (n_cells + 1) && !done; i++){
456 >      for (j = 0; j < (n_cells + 1) && !done; j++){
457 >        if (i < n_cells){
458 >          if (j < n_cells){
459 >            start_ndx = n_cells;
460 >          }
461 >          else
462 >            start_ndx = 0;
463 >        }
464 >        else
465 >          start_ndx = 0;
466  
467 <  if( i < n_cells ){
467 >        for (k = start_ndx; k < (n_cells + 1) && !done; k++){
468 >          makeElement(i * cellx, j * celly, k * cellz);
469 >          done = (current_mol >= tot_nmol);
470  
471 <    if( j < n_cells ){
472 <      start_ndx = n_cells;
473 <    }
474 <    else start_ndx = 0;
475 <  }
495 <  else start_ndx = 0;
471 >          if (!done && n_per_extra > 1){
472 >            makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly,
473 >                        k * cellz);
474 >            done = (current_mol >= tot_nmol);
475 >          }
476  
477 <  for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
478 <
479 <    makeElement( i * cellx,
480 <           j * celly,
481 <           k * cellz );
502 <    done = ( current_mol >= tot_nmol );
503 <
504 <    if( !done && n_per_extra > 1 ){
505 <      makeElement( i * cellx + 0.5 * cellx,
506 <       j * celly + 0.5 * celly,
507 <       k * cellz );
508 <      done = ( current_mol >= tot_nmol );
509 <    }
510 <
511 <    if( !done && n_per_extra > 2){
512 <      makeElement( i * cellx,
513 <       j * celly + 0.5 * celly,
514 <       k * cellz + 0.5 * cellz );
515 <      done = ( current_mol >= tot_nmol );
516 <    }
477 >          if (!done && n_per_extra > 2){
478 >            makeElement(i * cellx, j * celly + 0.5 * celly,
479 >                        k * cellz + 0.5 * cellz);
480 >            done = (current_mol >= tot_nmol);
481 >          }
482  
483 <    if( !done && n_per_extra > 3){
484 <      makeElement( i * cellx + 0.5 * cellx,
485 <       j * celly,
486 <       k * cellz + 0.5 * cellz );
487 <      done = ( current_mol >= tot_nmol );
488 <    }
524 <  }
483 >          if (!done && n_per_extra > 3){
484 >            makeElement(i * cellx + 0.5 * cellx, j * celly,
485 >                        k * cellz + 0.5 * cellz);
486 >            done = (current_mol >= tot_nmol);
487 >          }
488 >        }
489        }
490      }
491    }
492  
493 <  for( i=0; i<info[0].n_atoms; i++ ){
494 <    info[0].atoms[i]->setVel( vel );
493 >  for (i = 0; i < info[0].n_atoms; i++){
494 >    info[0].atoms[i]->setVel(vel);
495    }
496   }
497  
498 < void SimSetup::makeElement( double x, double y, double z ){
535 <
498 > void SimSetup::makeElement(double x, double y, double z){
499    int k;
500    AtomStamp* current_atom;
501    DirectionalAtom* dAtom;
502    double rotMat[3][3];
503    double pos[3];
504  
505 <  for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
506 <
507 <    current_atom = comp_stamps[current_comp]->getAtom( k );
508 <    if( !current_atom->havePosition() ){
509 <      sprintf( painCave.errMsg,
510 <         "SimSetup:initFromBass error.\n"
511 <         "\tComponent %s, atom %s does not have a position specified.\n"
512 <         "\tThe initialization routine is unable to give a start"
513 <         " position.\n",
551 <         comp_stamps[current_comp]->getID(),
552 <         current_atom->getType() );
505 >  for (k = 0; k < comp_stamps[current_comp]->getNAtoms(); k++){
506 >    current_atom = comp_stamps[current_comp]->getAtom(k);
507 >    if (!current_atom->havePosition()){
508 >      sprintf(painCave.errMsg,
509 >              "SimSetup:initFromBass error.\n"
510 >              "\tComponent %s, atom %s does not have a position specified.\n"
511 >              "\tThe initialization routine is unable to give a start"
512 >              " position.\n",
513 >              comp_stamps[current_comp]->getID(), current_atom->getType());
514        painCave.isFatal = 1;
515        simError();
516      }
517 <    
517 >
518      pos[0] = x + current_atom->getPosX();
519      pos[1] = y + current_atom->getPosY();
520      pos[2] = z + current_atom->getPosZ();
560    
561    info[0].atoms[current_atom_ndx]->setPos( pos );
521  
522 <    if( info[0].atoms[current_atom_ndx]->isDirectional() ){
522 >    info[0].atoms[current_atom_ndx]->setPos(pos);
523  
524 <      dAtom = (DirectionalAtom *)info[0].atoms[current_atom_ndx];
524 >    if (info[0].atoms[current_atom_ndx]->isDirectional()){
525 >      dAtom = (DirectionalAtom *) info[0].atoms[current_atom_ndx];
526  
527        rotMat[0][0] = 1.0;
528        rotMat[0][1] = 0.0;
# Line 576 | Line 536 | void SimSetup::makeElement( double x, double y, double
536        rotMat[2][1] = 0.0;
537        rotMat[2][2] = 1.0;
538  
539 <      dAtom->setA( rotMat );
539 >      dAtom->setA(rotMat);
540      }
541  
542      current_atom_ndx++;
# Line 585 | Line 545 | void SimSetup::makeElement( double x, double y, double
545    current_mol++;
546    current_comp_mol++;
547  
548 <  if( current_comp_mol >= components_nmol[current_comp] ){
589 <
548 >  if (current_comp_mol >= components_nmol[current_comp]){
549      current_comp_mol = 0;
550      current_comp++;
551    }
552   }
553  
554  
555 < void SimSetup::gatherInfo( void ){
556 <  int i,j,k;
555 > void SimSetup::gatherInfo(void){
556 >  int i, j, k;
557  
558    ensembleCase = -1;
559    ffCase = -1;
560  
561    // set the easy ones first
562  
563 <  for( i=0; i<nInfo; i++){
563 >  for (i = 0; i < nInfo; i++){
564      info[i].target_temp = globals->getTargetTemp();
565      info[i].dt = globals->getDt();
566      info[i].run_time = globals->getRunTime();
# Line 611 | Line 570 | void SimSetup::gatherInfo( void ){
570  
571    // get the forceField
572  
573 <  strcpy( force_field, globals->getForceField() );
573 >  strcpy(force_field, globals->getForceField());
574  
575 <  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
576 <  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
577 <  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
575 >  if (!strcasecmp(force_field, "DUFF")){
576 >    ffCase = FF_DUFF;
577 >  }
578 >  else if (!strcasecmp(force_field, "LJ")){
579 >    ffCase = FF_LJ;
580 >  }
581 >  else if (!strcasecmp(force_field, "EAM")){
582 >    ffCase = FF_EAM;
583 >  }
584    else{
585 <    sprintf( painCave.errMsg,
586 <       "SimSetup Error. Unrecognized force field -> %s\n",
587 <       force_field );
588 <    painCave.isFatal = 1;
624 <    simError();
585 >    sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
586 >            force_field);
587 >         painCave.isFatal = 1;
588 >         simError();
589    }
590  
591 <  // get the ensemble
591 >    // get the ensemble
592  
593 <  strcpy( ensemble, globals->getEnsemble() );
593 >  strcpy(ensemble, globals->getEnsemble());
594  
595 <  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
596 <  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
597 <  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )
595 >  if (!strcasecmp(ensemble, "NVE")){
596 >    ensembleCase = NVE_ENS;
597 >  }
598 >  else if (!strcasecmp(ensemble, "NVT")){
599 >    ensembleCase = NVT_ENS;
600 >  }
601 >  else if (!strcasecmp(ensemble, "NPTi") || !strcasecmp(ensemble, "NPT")){
602      ensembleCase = NPTi_ENS;
603 <  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
604 <  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
605 <  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
603 >  }
604 >  else if (!strcasecmp(ensemble, "NPTf")){
605 >    ensembleCase = NPTf_ENS;
606 >  }
607 >  else if (!strcasecmp(ensemble, "NPTim")){
608 >    ensembleCase = NPTim_ENS;
609 >  }
610 >  else if (!strcasecmp(ensemble, "NPTfm")){
611 >    ensembleCase = NPTfm_ENS;
612 >  }
613    else{
614 <    sprintf( painCave.errMsg,
615 <       "SimSetup Warning. Unrecognized Ensemble -> %s, "
616 <             "reverting to NVE for this simulation.\n",
617 <       ensemble );
618 <    painCave.isFatal = 0;
619 <    simError();
620 <    strcpy( ensemble, "NVE" );
621 <    ensembleCase = NVE_ENS;
614 >    sprintf(painCave.errMsg,
615 >            "SimSetup Warning. Unrecognized Ensemble -> %s, "
616 >            "reverting to NVE for this simulation.\n",
617 >            ensemble);
618 >         painCave.isFatal = 0;
619 >         simError();
620 >         strcpy(ensemble, "NVE");
621 >         ensembleCase = NVE_ENS;
622    }  
648  
649  for(i=0; i<nInfo; i++){
650    
651    strcpy( info[i].ensemble, ensemble );
623  
624 +  for (i = 0; i < nInfo; i++){
625 +    strcpy(info[i].ensemble, ensemble);
626 +
627      // get the mixing rule
628  
629 <    strcpy( info[i].mixingRule, globals->getMixingRule() );
629 >    strcpy(info[i].mixingRule, globals->getMixingRule());
630      info[i].usePBC = globals->getPBC();
631    }
632 <  
632 >
633    // get the components and calculate the tot_nMol and indvidual n_mol
634 <
634 >
635    the_components = globals->getComponents();
636    components_nmol = new int[n_components];
637  
638  
639 <  if( !globals->haveNMol() ){
639 >  if (!globals->haveNMol()){
640      // we don't have the total number of molecules, so we assume it is
641      // given in each component
642  
643      tot_nmol = 0;
644 <    for( i=0; i<n_components; i++ ){
645 <
646 <      if( !the_components[i]->haveNMol() ){
647 <  // we have a problem
648 <  sprintf( painCave.errMsg,
649 <     "SimSetup Error. No global NMol or component NMol"
650 <     " given. Cannot calculate the number of atoms.\n" );
651 <  painCave.isFatal = 1;
678 <  simError();
644 >    for (i = 0; i < n_components; i++){
645 >      if (!the_components[i]->haveNMol()){
646 >        // we have a problem
647 >        sprintf(painCave.errMsg,
648 >                "SimSetup Error. No global NMol or component NMol"
649 >                " given. Cannot calculate the number of atoms.\n");
650 >        painCave.isFatal = 1;
651 >        simError();
652        }
653  
654        tot_nmol += the_components[i]->getNMol();
# Line 683 | Line 656 | void SimSetup::gatherInfo( void ){
656      }
657    }
658    else{
659 <    sprintf( painCave.errMsg,
660 <       "SimSetup error.\n"
661 <       "\tSorry, the ability to specify total"
662 <       " nMols and then give molfractions in the components\n"
663 <       "\tis not currently supported."
664 <       " Please give nMol in the components.\n" );
659 >    sprintf(painCave.errMsg,
660 >            "SimSetup error.\n"
661 >            "\tSorry, the ability to specify total"
662 >            " nMols and then give molfractions in the components\n"
663 >            "\tis not currently supported."
664 >            " Please give nMol in the components.\n");
665      painCave.isFatal = 1;
666      simError();
667    }
668  
669    // set the status, sample, and thermal kick times
697  
698  for(i=0; i<nInfo; i++){
670  
671 <    if( globals->haveSampleTime() ){
671 >  for (i = 0; i < nInfo; i++){
672 >    if (globals->haveSampleTime()){
673        info[i].sampleTime = globals->getSampleTime();
674        info[i].statusTime = info[i].sampleTime;
675        info[i].thermalTime = info[i].sampleTime;
# Line 707 | Line 679 | void SimSetup::gatherInfo( void ){
679        info[i].statusTime = info[i].sampleTime;
680        info[i].thermalTime = info[i].sampleTime;
681      }
682 <    
683 <    if( globals->haveStatusTime() ){
682 >
683 >    if (globals->haveStatusTime()){
684        info[i].statusTime = globals->getStatusTime();
685      }
686 <    
687 <    if( globals->haveThermalTime() ){
686 >
687 >    if (globals->haveThermalTime()){
688        info[i].thermalTime = globals->getThermalTime();
689      }
690  
691      // check for the temperature set flag
692  
693 <    if( globals->haveTempSet() ) info[i].setTemp = globals->getTempSet();
694 <    
693 >    if (globals->haveTempSet())
694 >      info[i].setTemp = globals->getTempSet();
695 >
696      // get some of the tricky things that may still be in the globals
697 <    
697 >
698      double boxVector[3];
699 <    if( globals->haveBox() ){
699 >    if (globals->haveBox()){
700        boxVector[0] = globals->getBox();
701        boxVector[1] = globals->getBox();
702        boxVector[2] = globals->getBox();
703 <      
704 <      info[i].setBox( boxVector );
703 >
704 >      info[i].setBox(boxVector);
705      }
706 <    else if( globals->haveDensity() ){
734 <      
706 >    else if (globals->haveDensity()){
707        double vol;
708 <      vol = (double)tot_nmol / globals->getDensity();
709 <      boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
708 >      vol = (double) tot_nmol / globals->getDensity();
709 >      boxVector[0] = pow(vol, (1.0 / 3.0));
710        boxVector[1] = boxVector[0];
711        boxVector[2] = boxVector[0];
712 <      
713 <      info[i].setBox( boxVector );
714 <  }
712 >
713 >      info[i].setBox(boxVector);
714 >    }
715      else{
716 <      if( !globals->haveBoxX() ){
717 <  sprintf( painCave.errMsg,
718 <     "SimSetup error, no periodic BoxX size given.\n" );
719 <  painCave.isFatal = 1;
720 <  simError();
716 >      if (!globals->haveBoxX()){
717 >        sprintf(painCave.errMsg,
718 >                "SimSetup error, no periodic BoxX size given.\n");
719 >        painCave.isFatal = 1;
720 >        simError();
721        }
722        boxVector[0] = globals->getBoxX();
723 <      
724 <      if( !globals->haveBoxY() ){
725 <  sprintf( painCave.errMsg,
726 <     "SimSetup error, no periodic BoxY size given.\n" );
727 <  painCave.isFatal = 1;
728 <  simError();
723 >
724 >      if (!globals->haveBoxY()){
725 >        sprintf(painCave.errMsg,
726 >                "SimSetup error, no periodic BoxY size given.\n");
727 >        painCave.isFatal = 1;
728 >        simError();
729        }
730        boxVector[1] = globals->getBoxY();
731 <      
732 <      if( !globals->haveBoxZ() ){
733 <  sprintf( painCave.errMsg,
734 <     "SimSetup error, no periodic BoxZ size given.\n" );
735 <  painCave.isFatal = 1;
736 <  simError();
731 >
732 >      if (!globals->haveBoxZ()){
733 >        sprintf(painCave.errMsg,
734 >                "SimSetup error, no periodic BoxZ size given.\n");
735 >        painCave.isFatal = 1;
736 >        simError();
737        }
738        boxVector[2] = globals->getBoxZ();
739 <      
740 <      info[i].setBox( boxVector );
739 >
740 >      info[i].setBox(boxVector);
741      }
742    }
743  
744 +  //setup seed for random number generator
745    int seedValue;
773  struct timeval now_time_val;
774  struct timezone time_zone;
746  
747 <  if(globals->haveSeed()){
747 >  if (globals->haveSeed()){
748      seedValue = globals->getSeed();
749 <  }
749 >
750 >    if(seedValue / 1E9 == 0){
751 >      sprintf(painCave.errMsg,
752 >              "Seed for sprng library should contain at least 9 digits\n"
753 >              "OOPSE will generate a seed for user\n");
754 >      painCave.isFatal = 0;
755 >      simError();
756 >
757 >      //using seed generated by system instead of invalid seed set by user
758 > #ifndef IS_MPI
759 >      seedValue = make_sprng_seed();
760 > #else
761 >      if (worldRank == 0){
762 >        seedValue = make_sprng_seed();
763 >      }
764 >      MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
765 > #endif      
766 >    }
767 >  }//end of if branch of globals->haveSeed()
768    else{
769 +    
770   #ifndef IS_MPI
771 <    gettimeofday(&now_time_val, &time_zone);  // get the time now  
782 <    seedValue = (int) now_time_val.tv_sec;      //  convert to epoch time
771 >    seedValue = make_sprng_seed();
772   #else
773 <    if(worldRank == 0){
774 <      gettimeofday(&now_time_val, &time_zone);  // get the time now
786 <      seedValue = (int) now_time_val.tv_sec;      //  convert to epoch time      
773 >    if (worldRank == 0){
774 >      seedValue = make_sprng_seed();
775      }
776 <     MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
776 >    MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
777   #endif
778 <  }
778 >  }//end of globals->haveSeed()
779  
780 <  for(int i = 0; i < nInfo; i++){
780 >  for (int i = 0; i < nInfo; i++){
781      info[i].setSeed(seedValue);
782    }
783  
784   #ifdef IS_MPI
785 <  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
785 >  strcpy(checkPointMsg, "Succesfully gathered all information from Bass\n");
786    MPIcheckPoint();
787   #endif // is_mpi
800
788   }
789  
790  
791 < void SimSetup::finalInfoCheck( void ){
791 > void SimSetup::finalInfoCheck(void){
792    int index;
793    int usesDipoles;
794    int i;
795  
796 <  for(i=0; i<nInfo; i++){
796 >  for (i = 0; i < nInfo; i++){
797      // check electrostatic parameters
798 <    
798 >
799      index = 0;
800      usesDipoles = 0;
801 <    while( (index < info[i].n_atoms) && !usesDipoles ){
801 >    while ((index < info[i].n_atoms) && !usesDipoles){
802        usesDipoles = (info[i].atoms[index])->hasDipole();
803        index++;
804      }
805 <    
805 >
806   #ifdef IS_MPI
807      int myUse = usesDipoles;
808 <    MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
808 >    MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
809   #endif //is_mpi
810 <    
810 >
811      double theEcr, theEst;
812 <  
813 <    if (globals->getUseRF() ) {
812 >
813 >    if (globals->getUseRF()){
814        info[i].useReactionField = 1;
815 <      
816 <      if( !globals->haveECR() ){
817 <  sprintf( painCave.errMsg,
818 <     "SimSetup Warning: using default value of 1/2 the smallest "
819 <     "box length for the electrostaticCutoffRadius.\n"
820 <     "I hope you have a very fast processor!\n");
821 <  painCave.isFatal = 0;
822 <  simError();
823 <  double smallest;
824 <  smallest = info[i].boxL[0];
825 <  if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
826 <  if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
827 <  theEcr = 0.5 * smallest;
828 <      } else {
829 <  theEcr = globals->getECR();
815 >
816 >      if (!globals->haveECR()){
817 >        sprintf(painCave.errMsg,
818 >                "SimSetup Warning: using default value of 1/2 the smallest "
819 >                "box length for the electrostaticCutoffRadius.\n"
820 >                "I hope you have a very fast processor!\n");
821 >        painCave.isFatal = 0;
822 >        simError();
823 >        double smallest;
824 >        smallest = info[i].boxL[0];
825 >        if (info[i].boxL[1] <= smallest)
826 >          smallest = info[i].boxL[1];
827 >        if (info[i].boxL[2] <= smallest)
828 >          smallest = info[i].boxL[2];
829 >        theEcr = 0.5 * smallest;
830        }
831 <      
832 <      if( !globals->haveEST() ){
846 <  sprintf( painCave.errMsg,
847 <     "SimSetup Warning: using default value of 0.05 * the "
848 <     "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
849 <     );
850 <  painCave.isFatal = 0;
851 <  simError();
852 <  theEst = 0.05 * theEcr;
853 <      } else {
854 <  theEst= globals->getEST();
831 >      else{
832 >        theEcr = globals->getECR();
833        }
834 <      
835 <      info[i].setEcr( theEcr, theEst );
836 <      
837 <      if(!globals->haveDielectric() ){
838 <  sprintf( painCave.errMsg,
839 <     "SimSetup Error: You are trying to use Reaction Field without"
840 <     "setting a dielectric constant!\n"
841 <     );
864 <  painCave.isFatal = 1;
865 <  simError();
834 >
835 >      if (!globals->haveEST()){
836 >        sprintf(painCave.errMsg,
837 >                "SimSetup Warning: using default value of 0.05 * the "
838 >                "electrostaticCutoffRadius for the electrostaticSkinThickness\n");
839 >        painCave.isFatal = 0;
840 >        simError();
841 >        theEst = 0.05 * theEcr;
842        }
843 <      info[i].dielectric = globals->getDielectric();  
844 <    }
869 <    else {
870 <      if (usesDipoles) {
871 <  
872 <  if( !globals->haveECR() ){
873 <    sprintf( painCave.errMsg,
874 <       "SimSetup Warning: using default value of 1/2 the smallest "
875 <       "box length for the electrostaticCutoffRadius.\n"
876 <       "I hope you have a very fast processor!\n");
877 <    painCave.isFatal = 0;
878 <    simError();
879 <    double smallest;
880 <    smallest = info[i].boxL[0];
881 <    if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
882 <    if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
883 <    theEcr = 0.5 * smallest;
884 <  } else {
885 <    theEcr = globals->getECR();
886 <  }
887 <  
888 <  if( !globals->haveEST() ){
889 <    sprintf( painCave.errMsg,
890 <       "SimSetup Warning: using default value of 0.05 * the "
891 <       "electrostaticCutoffRadius for the "
892 <       "electrostaticSkinThickness\n"
893 <       );
894 <    painCave.isFatal = 0;
895 <    simError();
896 <    theEst = 0.05 * theEcr;
897 <  } else {
898 <    theEst= globals->getEST();
899 <  }
900 <  
901 <  info[i].setEcr( theEcr, theEst );
843 >      else{
844 >        theEst = globals->getEST();
845        }
846 <    }  
846 >
847 >      info[i].setEcr(theEcr, theEst);
848 >
849 >      if (!globals->haveDielectric()){
850 >        sprintf(painCave.errMsg,
851 >                "SimSetup Error: You are trying to use Reaction Field without"
852 >                "setting a dielectric constant!\n");
853 >        painCave.isFatal = 1;
854 >        simError();
855 >      }
856 >      info[i].dielectric = globals->getDielectric();
857 >    }
858 >    else{
859 >      if (usesDipoles){
860 >        if (!globals->haveECR()){
861 >          sprintf(painCave.errMsg,
862 >                  "SimSetup Warning: using default value of 1/2 the smallest "
863 >                  "box length for the electrostaticCutoffRadius.\n"
864 >                  "I hope you have a very fast processor!\n");
865 >          painCave.isFatal = 0;
866 >          simError();
867 >          double smallest;
868 >          smallest = info[i].boxL[0];
869 >          if (info[i].boxL[1] <= smallest)
870 >            smallest = info[i].boxL[1];
871 >          if (info[i].boxL[2] <= smallest)
872 >            smallest = info[i].boxL[2];
873 >          theEcr = 0.5 * smallest;
874 >        }
875 >        else{
876 >          theEcr = globals->getECR();
877 >        }
878 >
879 >        if (!globals->haveEST()){
880 >          sprintf(painCave.errMsg,
881 >                  "SimSetup Warning: using default value of 0.05 * the "
882 >                  "electrostaticCutoffRadius for the "
883 >                  "electrostaticSkinThickness\n");
884 >          painCave.isFatal = 0;
885 >          simError();
886 >          theEst = 0.05 * theEcr;
887 >        }
888 >        else{
889 >          theEst = globals->getEST();
890 >        }
891 >
892 >        info[i].setEcr(theEcr, theEst);
893 >      }
894 >    }
895    }
896  
897   #ifdef IS_MPI
898 <  strcpy( checkPointMsg, "post processing checks out" );
898 >  strcpy(checkPointMsg, "post processing checks out");
899    MPIcheckPoint();
900   #endif // is_mpi
910
901   }
902  
903 < void SimSetup::initSystemCoords( void ){
903 > void SimSetup::initSystemCoords(void){
904    int i;
905 <  
905 >
906    char* inName;
907  
908 <  (info[0].getConfiguration())->createArrays( info[0].n_atoms );
908 >  (info[0].getConfiguration())->createArrays(info[0].n_atoms);
909  
910 <  for(i=0; i<info[0].n_atoms; i++) info[0].atoms[i]->setCoords();
911 <  
912 <  if( globals->haveInitialConfig() ){
913 <    
914 <    InitializeFromFile* fileInit;
915 < #ifdef IS_MPI // is_mpi
916 <    if( worldRank == 0 ){
910 >  for (i = 0; i < info[0].n_atoms; i++)
911 >    info[0].atoms[i]->setCoords();
912 >
913 >  if (globals->haveInitialConfig()){
914 >    InitializeFromFile* fileInit;
915 > #ifdef IS_MPI // is_mpi
916 >    if (worldRank == 0){
917   #endif //is_mpi
918        inName = globals->getInitialConfig();
919        double* tempDouble = new double[1000000];
920 <      fileInit = new InitializeFromFile( inName );
920 >      fileInit = new InitializeFromFile(inName);
921   #ifdef IS_MPI
922 <    }else fileInit = new InitializeFromFile( NULL );
922 >    }
923 >    else
924 >      fileInit = new InitializeFromFile(NULL);
925   #endif
926 <    fileInit->readInit( info ); // default velocities on
927 <    
926 >    fileInit->readInit(info); // default velocities on
927 >
928      delete fileInit;
929    }
930    else{
939    
931   #ifdef IS_MPI
932 <    
932 >
933      // no init from bass
934 <    
935 <    sprintf( painCave.errMsg,
936 <       "Cannot intialize a parallel simulation without an initial configuration file.\n" );
934 >
935 >    sprintf(painCave.errMsg,
936 >            "Cannot intialize a parallel simulation without an initial configuration file.\n");
937      painCave.isFatal;
938      simError();
939 <    
939 >
940   #else
941 <    
941 >
942      initFromBass();
943 <    
944 <    
943 >
944 >
945   #endif
946    }
947 <  
947 >
948   #ifdef IS_MPI
949 <  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
949 >  strcpy(checkPointMsg, "Successfully read in the initial configuration");
950    MPIcheckPoint();
951   #endif // is_mpi
961  
952   }
953  
954  
955 < void SimSetup::makeOutNames( void ){
966 <  
955 > void SimSetup::makeOutNames(void){
956    int k;
957  
969  
970  for(k=0; k<nInfo; k++){
958  
959 +  for (k = 0; k < nInfo; k++){
960   #ifdef IS_MPI
961 <    if( worldRank == 0 ){
961 >    if (worldRank == 0){
962   #endif // is_mpi
963 <      
964 <      if( globals->haveFinalConfig() ){
965 <  strcpy( info[k].finalName, globals->getFinalConfig() );
963 >
964 >      if (globals->haveFinalConfig()){
965 >        strcpy(info[k].finalName, globals->getFinalConfig());
966        }
967        else{
968 <  strcpy( info[k].finalName, inFileName );
969 <  char* endTest;
970 <  int nameLength = strlen( info[k].finalName );
971 <  endTest = &(info[k].finalName[nameLength - 5]);
972 <  if( !strcmp( endTest, ".bass" ) ){
973 <    strcpy( endTest, ".eor" );
974 <  }
975 <  else if( !strcmp( endTest, ".BASS" ) ){
976 <    strcpy( endTest, ".eor" );
977 <  }
978 <  else{
979 <    endTest = &(info[k].finalName[nameLength - 4]);
980 <    if( !strcmp( endTest, ".bss" ) ){
981 <      strcpy( endTest, ".eor" );
982 <    }
983 <    else if( !strcmp( endTest, ".mdl" ) ){
984 <      strcpy( endTest, ".eor" );
985 <    }
986 <    else{
987 <      strcat( info[k].finalName, ".eor" );
988 <    }
989 <  }
968 >        strcpy(info[k].finalName, inFileName);
969 >        char* endTest;
970 >        int nameLength = strlen(info[k].finalName);
971 >        endTest = &(info[k].finalName[nameLength - 5]);
972 >        if (!strcmp(endTest, ".bass")){
973 >          strcpy(endTest, ".eor");
974 >        }
975 >        else if (!strcmp(endTest, ".BASS")){
976 >          strcpy(endTest, ".eor");
977 >        }
978 >        else{
979 >          endTest = &(info[k].finalName[nameLength - 4]);
980 >          if (!strcmp(endTest, ".bss")){
981 >            strcpy(endTest, ".eor");
982 >          }
983 >          else if (!strcmp(endTest, ".mdl")){
984 >            strcpy(endTest, ".eor");
985 >          }
986 >          else{
987 >            strcat(info[k].finalName, ".eor");
988 >          }
989 >        }
990        }
991 <      
991 >
992        // make the sample and status out names
993 <      
994 <      strcpy( info[k].sampleName, inFileName );
993 >
994 >      strcpy(info[k].sampleName, inFileName);
995        char* endTest;
996 <      int nameLength = strlen( info[k].sampleName );
996 >      int nameLength = strlen(info[k].sampleName);
997        endTest = &(info[k].sampleName[nameLength - 5]);
998 <      if( !strcmp( endTest, ".bass" ) ){
999 <  strcpy( endTest, ".dump" );
998 >      if (!strcmp(endTest, ".bass")){
999 >        strcpy(endTest, ".dump");
1000        }
1001 <      else if( !strcmp( endTest, ".BASS" ) ){
1002 <  strcpy( endTest, ".dump" );
1001 >      else if (!strcmp(endTest, ".BASS")){
1002 >        strcpy(endTest, ".dump");
1003        }
1004        else{
1005 <  endTest = &(info[k].sampleName[nameLength - 4]);
1006 <  if( !strcmp( endTest, ".bss" ) ){
1007 <    strcpy( endTest, ".dump" );
1008 <  }
1009 <  else if( !strcmp( endTest, ".mdl" ) ){
1010 <    strcpy( endTest, ".dump" );
1011 <  }
1012 <  else{
1013 <    strcat( info[k].sampleName, ".dump" );
1014 <  }
1005 >        endTest = &(info[k].sampleName[nameLength - 4]);
1006 >        if (!strcmp(endTest, ".bss")){
1007 >          strcpy(endTest, ".dump");
1008 >        }
1009 >        else if (!strcmp(endTest, ".mdl")){
1010 >          strcpy(endTest, ".dump");
1011 >        }
1012 >        else{
1013 >          strcat(info[k].sampleName, ".dump");
1014 >        }
1015        }
1016 <      
1017 <      strcpy( info[k].statusName, inFileName );
1018 <      nameLength = strlen( info[k].statusName );
1016 >
1017 >      strcpy(info[k].statusName, inFileName);
1018 >      nameLength = strlen(info[k].statusName);
1019        endTest = &(info[k].statusName[nameLength - 5]);
1020 <      if( !strcmp( endTest, ".bass" ) ){
1021 <  strcpy( endTest, ".stat" );
1020 >      if (!strcmp(endTest, ".bass")){
1021 >        strcpy(endTest, ".stat");
1022        }
1023 <      else if( !strcmp( endTest, ".BASS" ) ){
1024 <  strcpy( endTest, ".stat" );
1023 >      else if (!strcmp(endTest, ".BASS")){
1024 >        strcpy(endTest, ".stat");
1025        }
1026        else{
1027 <  endTest = &(info[k].statusName[nameLength - 4]);
1028 <  if( !strcmp( endTest, ".bss" ) ){
1029 <    strcpy( endTest, ".stat" );
1030 <  }
1031 <  else if( !strcmp( endTest, ".mdl" ) ){
1032 <    strcpy( endTest, ".stat" );
1033 <  }
1034 <  else{
1035 <    strcat( info[k].statusName, ".stat" );
1036 <  }
1027 >        endTest = &(info[k].statusName[nameLength - 4]);
1028 >        if (!strcmp(endTest, ".bss")){
1029 >          strcpy(endTest, ".stat");
1030 >        }
1031 >        else if (!strcmp(endTest, ".mdl")){
1032 >          strcpy(endTest, ".stat");
1033 >        }
1034 >        else{
1035 >          strcat(info[k].statusName, ".stat");
1036 >        }
1037        }
1038 <      
1038 >
1039   #ifdef IS_MPI
1040 +
1041      }
1042   #endif // is_mpi
1043    }
1044   }
1045  
1046  
1047 < void SimSetup::sysObjectsCreation( void ){
1048 <  
1049 <  int i,k;
1061 <  
1047 > void SimSetup::sysObjectsCreation(void){
1048 >  int i, k;
1049 >
1050    // create the forceField
1051  
1052    createFF();
# Line 1073 | Line 1061 | void SimSetup::sysObjectsCreation( void ){
1061  
1062   #ifdef IS_MPI
1063    // divide the molecules among the processors
1064 <  
1064 >
1065    mpiMolDivide();
1066   #endif //is_mpi
1067 <  
1067 >
1068    // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
1069 <
1069 >
1070    makeSysArrays();
1071  
1072    // make and initialize the molecules (all but atomic coordinates)
1073 <
1073 >
1074    makeMolecules();
1075 <  
1076 <  for(k=0; k<nInfo; k++){
1075 >
1076 >  for (k = 0; k < nInfo; k++){
1077      info[k].identArray = new int[info[k].n_atoms];
1078 <    for(i=0; i<info[k].n_atoms; i++){
1078 >    for (i = 0; i < info[k].n_atoms; i++){
1079        info[k].identArray[i] = info[k].atoms[i]->getIdent();
1080      }
1081    }
1082   }
1083  
1084  
1085 < void SimSetup::createFF( void ){
1085 > void SimSetup::createFF(void){
1086 >  switch (ffCase){
1087 >    case FF_DUFF:
1088 >      the_ff = new DUFF();
1089 >      break;
1090  
1091 <  switch( ffCase ){
1091 >    case FF_LJ:
1092 >      the_ff = new LJFF();
1093 >      break;
1094  
1095 <  case FF_DUFF:
1096 <    the_ff = new DUFF();
1097 <    break;
1095 >    case FF_EAM:
1096 >      the_ff = new EAM_FF();
1097 >      break;
1098  
1099 <  case FF_LJ:
1100 <    the_ff = new LJFF();
1101 <    break;
1102 <
1103 <  case FF_EAM:
1110 <    the_ff = new EAM_FF();
1111 <    break;
1112 <
1113 <  default:
1114 <    sprintf( painCave.errMsg,
1115 <       "SimSetup Error. Unrecognized force field in case statement.\n");
1116 <    painCave.isFatal = 1;
1117 <    simError();
1099 >    default:
1100 >      sprintf(painCave.errMsg,
1101 >              "SimSetup Error. Unrecognized force field in case statement.\n");
1102 >      painCave.isFatal = 1;
1103 >      simError();
1104    }
1105  
1106   #ifdef IS_MPI
1107 <  strcpy( checkPointMsg, "ForceField creation successful" );
1107 >  strcpy(checkPointMsg, "ForceField creation successful");
1108    MPIcheckPoint();
1109   #endif // is_mpi
1124
1110   }
1111  
1112  
1113 < void SimSetup::compList( void ){
1129 <
1113 > void SimSetup::compList(void){
1114    int i;
1115    char* id;
1116    LinkedMolStamp* headStamp = new LinkedMolStamp();
1117    LinkedMolStamp* currentStamp = NULL;
1118 <  comp_stamps = new MoleculeStamp*[n_components];
1119 <  
1118 >  comp_stamps = new MoleculeStamp * [n_components];
1119 >
1120    // make an array of molecule stamps that match the components used.
1121    // also extract the used stamps out into a separate linked list
1122 <  
1123 <  for(i=0; i<nInfo; i++){
1122 >
1123 >  for (i = 0; i < nInfo; i++){
1124      info[i].nComponents = n_components;
1125      info[i].componentsNmol = components_nmol;
1126      info[i].compStamps = comp_stamps;
1127      info[i].headStamp = headStamp;
1128    }
1145  
1129  
1147  for( i=0; i<n_components; i++ ){
1130  
1131 +  for (i = 0; i < n_components; i++){
1132      id = the_components[i]->getType();
1133      comp_stamps[i] = NULL;
1151    
1152    // check to make sure the component isn't already in the list
1134  
1135 <    comp_stamps[i] = headStamp->match( id );
1136 <    if( comp_stamps[i] == NULL ){
1137 <      
1135 >    // check to make sure the component isn't already in the list
1136 >
1137 >    comp_stamps[i] = headStamp->match(id);
1138 >    if (comp_stamps[i] == NULL){
1139        // extract the component from the list;
1140 <      
1141 <      currentStamp = stamps->extractMolStamp( id );
1142 <      if( currentStamp == NULL ){
1143 <  sprintf( painCave.errMsg,
1144 <     "SimSetup error: Component \"%s\" was not found in the "
1145 <     "list of declared molecules\n",
1146 <     id );
1147 <  painCave.isFatal = 1;
1148 <  simError();
1140 >
1141 >      currentStamp = stamps->extractMolStamp(id);
1142 >      if (currentStamp == NULL){
1143 >        sprintf(painCave.errMsg,
1144 >                "SimSetup error: Component \"%s\" was not found in the "
1145 >                "list of declared molecules\n",
1146 >                id);
1147 >        painCave.isFatal = 1;
1148 >        simError();
1149        }
1150 <      
1151 <      headStamp->add( currentStamp );
1152 <      comp_stamps[i] = headStamp->match( id );
1150 >
1151 >      headStamp->add(currentStamp);
1152 >      comp_stamps[i] = headStamp->match(id);
1153      }
1154    }
1155  
1156   #ifdef IS_MPI
1157 <  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
1157 >  strcpy(checkPointMsg, "Component stamps successfully extracted\n");
1158    MPIcheckPoint();
1159   #endif // is_mpi
1178
1179
1160   }
1161  
1162 < void SimSetup::calcSysValues( void ){
1162 > void SimSetup::calcSysValues(void){
1163    int i, j, k;
1164 <  
1165 <  int *molMembershipArray;
1166 <  
1164 >
1165 >  int* molMembershipArray;
1166 >
1167    tot_atoms = 0;
1168    tot_bonds = 0;
1169    tot_bends = 0;
1170    tot_torsions = 0;
1171 <  for( i=0; i<n_components; i++ ){
1172 <    
1173 <    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
1174 <    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
1195 <    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
1171 >  for (i = 0; i < n_components; i++){
1172 >    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1173 >    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1174 >    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1175      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1176    }
1177 <  
1177 >
1178    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1179    molMembershipArray = new int[tot_atoms];
1180 <  
1181 <  for(i=0; i<nInfo; i++){
1180 >
1181 >  for (i = 0; i < nInfo; i++){
1182      info[i].n_atoms = tot_atoms;
1183      info[i].n_bonds = tot_bonds;
1184      info[i].n_bends = tot_bends;
1185      info[i].n_torsions = tot_torsions;
1186      info[i].n_SRI = tot_SRI;
1187      info[i].n_mol = tot_nmol;
1188 <    
1188 >
1189      info[i].molMembershipArray = molMembershipArray;
1190 <  }
1190 >  }
1191   }
1192  
1193   #ifdef IS_MPI
1194  
1195 < void SimSetup::mpiMolDivide( void ){
1217 <  
1195 > void SimSetup::mpiMolDivide(void){
1196    int i, j, k;
1197    int localMol, allMol;
1198    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1199  
1200 <  mpiSim = new mpiSimulation( info );
1201 <  
1200 >  mpiSim = new mpiSimulation(info);
1201 >
1202    globalIndex = mpiSim->divideLabor();
1203  
1204    // set up the local variables
1205 <  
1205 >
1206    mol2proc = mpiSim->getMolToProcMap();
1207    molCompType = mpiSim->getMolComponentType();
1208 <  
1208 >
1209    allMol = 0;
1210    localMol = 0;
1211    local_atoms = 0;
# Line 1237 | Line 1215 | void SimSetup::mpiMolDivide( void ){
1215    globalAtomIndex = 0;
1216  
1217  
1218 <  for( i=0; i<n_components; i++ ){
1219 <
1220 <    for( j=0; j<components_nmol[i]; j++ ){
1221 <      
1222 <      if( mol2proc[allMol] == worldRank ){
1223 <  
1224 <  local_atoms +=    comp_stamps[i]->getNAtoms();
1225 <  local_bonds +=    comp_stamps[i]->getNBonds();
1248 <  local_bends +=    comp_stamps[i]->getNBends();
1249 <  local_torsions += comp_stamps[i]->getNTorsions();
1250 <  localMol++;
1218 >  for (i = 0; i < n_components; i++){
1219 >    for (j = 0; j < components_nmol[i]; j++){
1220 >      if (mol2proc[allMol] == worldRank){
1221 >        local_atoms += comp_stamps[i]->getNAtoms();
1222 >        local_bonds += comp_stamps[i]->getNBonds();
1223 >        local_bends += comp_stamps[i]->getNBends();
1224 >        local_torsions += comp_stamps[i]->getNTorsions();
1225 >        localMol++;
1226        }      
1227 <      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1227 >      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1228          info[0].molMembershipArray[globalAtomIndex] = allMol;
1229          globalAtomIndex++;
1230        }
1231  
1232 <      allMol++;      
1232 >      allMol++;
1233      }
1234    }
1235    local_SRI = local_bonds + local_bends + local_torsions;
1236 <  
1236 >
1237    info[0].n_atoms = mpiSim->getMyNlocal();  
1238 <  
1239 <  if( local_atoms != info[0].n_atoms ){
1240 <    sprintf( painCave.errMsg,
1241 <       "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1242 <       " localAtom (%d) are not equal.\n",
1243 <       info[0].n_atoms,
1269 <       local_atoms );
1238 >
1239 >  if (local_atoms != info[0].n_atoms){
1240 >    sprintf(painCave.errMsg,
1241 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1242 >            " localAtom (%d) are not equal.\n",
1243 >            info[0].n_atoms, local_atoms);
1244      painCave.isFatal = 1;
1245      simError();
1246    }
# Line 1277 | Line 1251 | void SimSetup::mpiMolDivide( void ){
1251    info[0].n_SRI = local_SRI;
1252    info[0].n_mol = localMol;
1253  
1254 <  strcpy( checkPointMsg, "Passed nlocal consistency check." );
1254 >  strcpy(checkPointMsg, "Passed nlocal consistency check.");
1255    MPIcheckPoint();
1256   }
1257 <
1257 >
1258   #endif // is_mpi
1259  
1260  
1261 < void SimSetup::makeSysArrays( void ){
1261 > void SimSetup::makeSysArrays(void){
1262    int i, j, k, l;
1263  
1264    Atom** the_atoms;
1265    Molecule* the_molecules;
1266    Exclude** the_excludes;
1267  
1268 <  
1269 <  for(l=0; l<nInfo; l++){
1296 <    
1268 >
1269 >  for (l = 0; l < nInfo; l++){
1270      // create the atom and short range interaction arrays
1271 <    
1272 <    the_atoms = new Atom*[info[l].n_atoms];
1271 >
1272 >    the_atoms = new Atom * [info[l].n_atoms];
1273      the_molecules = new Molecule[info[l].n_mol];
1274      int molIndex;
1275  
1276      // initialize the molecule's stampID's
1277 <    
1277 >
1278   #ifdef IS_MPI
1279 <    
1280 <    
1279 >
1280 >
1281      molIndex = 0;
1282 <    for(i=0; i<mpiSim->getTotNmol(); i++){
1283 <    
1284 <      if(mol2proc[i] == worldRank ){
1285 <  the_molecules[molIndex].setStampID( molCompType[i] );
1286 <  the_molecules[molIndex].setMyIndex( molIndex );
1287 <  the_molecules[molIndex].setGlobalIndex( i );
1315 <  molIndex++;
1282 >    for (i = 0; i < mpiSim->getTotNmol(); i++){
1283 >      if (mol2proc[i] == worldRank){
1284 >        the_molecules[molIndex].setStampID(molCompType[i]);
1285 >        the_molecules[molIndex].setMyIndex(molIndex);
1286 >        the_molecules[molIndex].setGlobalIndex(i);
1287 >        molIndex++;
1288        }
1289      }
1290 <    
1290 >
1291   #else // is_mpi
1292 <    
1292 >
1293      molIndex = 0;
1294      globalAtomIndex = 0;
1295 <    for(i=0; i<n_components; i++){
1296 <      for(j=0; j<components_nmol[i]; j++ ){
1297 <  the_molecules[molIndex].setStampID( i );
1298 <  the_molecules[molIndex].setMyIndex( molIndex );
1299 <  the_molecules[molIndex].setGlobalIndex( molIndex );
1300 <  for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1301 <    info[l].molMembershipArray[globalAtomIndex] = molIndex;
1302 <    globalAtomIndex++;
1303 <  }
1304 <  molIndex++;
1295 >    for (i = 0; i < n_components; i++){
1296 >      for (j = 0; j < components_nmol[i]; j++){
1297 >        the_molecules[molIndex].setStampID(i);
1298 >        the_molecules[molIndex].setMyIndex(molIndex);
1299 >        the_molecules[molIndex].setGlobalIndex(molIndex);
1300 >        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1301 >          info[l].molMembershipArray[globalAtomIndex] = molIndex;
1302 >          globalAtomIndex++;
1303 >        }
1304 >        molIndex++;
1305        }
1306      }
1307 <    
1308 <    
1307 >
1308 >
1309   #endif // is_mpi
1310  
1311  
1312 <    if( info[l].n_SRI ){
1341 <    
1312 >    if (info[l].n_SRI){
1313        Exclude::createArray(info[l].n_SRI);
1314 <      the_excludes = new Exclude*[info[l].n_SRI];
1315 <      for( int ex=0; ex<info[l].n_SRI; ex++){
1316 <  the_excludes[ex] = new Exclude(ex);
1314 >      the_excludes = new Exclude * [info[l].n_SRI];
1315 >      for (int ex = 0; ex < info[l].n_SRI; ex++){
1316 >        the_excludes[ex] = new Exclude(ex);
1317        }
1318        info[l].globalExcludes = new int;
1319        info[l].n_exclude = info[l].n_SRI;
1320      }
1321      else{
1322 <    
1323 <      Exclude::createArray( 1 );
1353 <      the_excludes = new Exclude*;
1322 >      Exclude::createArray(1);
1323 >      the_excludes = new Exclude * ;
1324        the_excludes[0] = new Exclude(0);
1325 <      the_excludes[0]->setPair( 0,0 );
1325 >      the_excludes[0]->setPair(0, 0);
1326        info[l].globalExcludes = new int;
1327        info[l].globalExcludes[0] = 0;
1328        info[l].n_exclude = 0;
# Line 1365 | Line 1335 | void SimSetup::makeSysArrays( void ){
1335      info[l].nGlobalExcludes = 0;
1336      info[l].excludes = the_excludes;
1337  
1338 <    the_ff->setSimInfo( info );
1369 <    
1338 >    the_ff->setSimInfo(info);
1339    }
1340   }
1341  
1342 < void SimSetup::makeIntegrator( void ){
1374 <
1342 > void SimSetup::makeIntegrator(void){
1343    int k;
1344  
1345 <  NVT<RealIntegrator>*  myNVT = NULL;
1345 >  NVT<RealIntegrator>* myNVT = NULL;
1346    NPTi<RealIntegrator>* myNPTi = NULL;
1347    NPTf<RealIntegrator>* myNPTf = NULL;
1348    NPTim<RealIntegrator>* myNPTim = NULL;
1349    NPTfm<RealIntegrator>* myNPTfm = NULL;
1382        
1383  for(k=0; k<nInfo; k++){
1384    
1385    switch( ensembleCase ){
1386      
1387    case NVE_ENS:
1388      if (globals->haveZconstraints()){
1389        setupZConstraint(info[k]);
1390        new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
1391     }
1350  
1351 <     else
1352 <        new NVE<RealIntegrator>( &(info[k]), the_ff );
1353 <      break;
1354 <      
1355 <    case NVT_ENS:
1356 <      if (globals->haveZconstraints()){
1357 <        setupZConstraint(info[k]);
1358 <        myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
1359 <      }
1360 <      else
1403 <        myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
1351 >  for (k = 0; k < nInfo; k++){
1352 >    switch (ensembleCase){
1353 >      case NVE_ENS:
1354 >        if (globals->haveZconstraints()){
1355 >          setupZConstraint(info[k]);
1356 >          new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1357 >        }
1358 >        else
1359 >          new NVE<RealIntegrator>(&(info[k]), the_ff);
1360 >        break;
1361  
1362 +      case NVT_ENS:
1363 +        if (globals->haveZconstraints()){
1364 +          setupZConstraint(info[k]);
1365 +          myNVT = new ZConstraint<NVT<RealIntegrator> >(&(info[k]), the_ff);
1366 +        }
1367 +        else
1368 +          myNVT = new NVT<RealIntegrator>(&(info[k]), the_ff);
1369 +
1370          myNVT->setTargetTemp(globals->getTargetTemp());
1371 <      
1372 <        if (globals->haveTauThermostat())
1371 >
1372 >        if (globals->haveTauThermostat())
1373            myNVT->setTauThermostat(globals->getTauThermostat());
1374 <      
1375 <        else {
1376 <          sprintf( painCave.errMsg,
1377 <                    "SimSetup error: If you use the NVT\n"
1413 <                    "    ensemble, you must set tauThermostat.\n");
1374 >        else{
1375 >          sprintf(painCave.errMsg,
1376 >                  "SimSetup error: If you use the NVT\n"
1377 >                  "    ensemble, you must set tauThermostat.\n");
1378            painCave.isFatal = 1;
1379            simError();
1380          }
1381          break;
1418      
1419    case NPTi_ENS:
1420      if (globals->haveZconstraints()){
1421             setupZConstraint(info[k]);
1422         myNPTi = new ZConstraint<NPTi<RealIntegrator> >( &(info[k]), the_ff );
1423      }
1424      else
1425        myNPTi = new NPTi<RealIntegrator>( &(info[k]), the_ff );
1382  
1383 <      myNPTi->setTargetTemp( globals->getTargetTemp() );
1384 <          
1385 <      if (globals->haveTargetPressure())
1386 <        myNPTi->setTargetPressure(globals->getTargetPressure());
1387 <      else {
1388 <         sprintf( painCave.errMsg,
1389 <                   "SimSetup error: If you use a constant pressure\n"
1434 <                   "    ensemble, you must set targetPressure in the BASS file.\n");
1435 <         painCave.isFatal = 1;
1436 <         simError();
1437 <      }
1438 <          
1439 <      if( globals->haveTauThermostat() )
1440 <        myNPTi->setTauThermostat( globals->getTauThermostat() );
1441 <      else{
1442 <         sprintf( painCave.errMsg,
1443 <                   "SimSetup error: If you use an NPT\n"
1444 <                  "    ensemble, you must set tauThermostat.\n");
1445 <         painCave.isFatal = 1;
1446 <         simError();
1447 <      }
1448 <          
1449 <      if( globals->haveTauBarostat() )
1450 <        myNPTi->setTauBarostat( globals->getTauBarostat() );
1451 <      else{
1452 <        sprintf( painCave.errMsg,
1453 <                  "SimSetup error: If you use an NPT\n"
1454 <                  "    ensemble, you must set tauBarostat.\n");
1455 <        painCave.isFatal = 1;
1456 <        simError();
1457 <       }
1458 <       break;
1459 <      
1460 <    case NPTf_ENS:
1461 <      if (globals->haveZconstraints()){
1462 <        setupZConstraint(info[k]);
1463 <        myNPTf = new ZConstraint<NPTf<RealIntegrator> >( &(info[k]), the_ff );
1464 <      }
1465 <      else
1466 <        myNPTf = new NPTf<RealIntegrator>( &(info[k]), the_ff );
1383 >      case NPTi_ENS:
1384 >        if (globals->haveZconstraints()){
1385 >          setupZConstraint(info[k]);
1386 >          myNPTi = new ZConstraint<NPTi<RealIntegrator> >(&(info[k]), the_ff);
1387 >        }
1388 >        else
1389 >          myNPTi = new NPTi<RealIntegrator>(&(info[k]), the_ff);
1390  
1391 <      myNPTf->setTargetTemp( globals->getTargetTemp());
1469 <          
1470 <      if (globals->haveTargetPressure())
1471 <        myNPTf->setTargetPressure(globals->getTargetPressure());
1472 <      else {
1473 <        sprintf( painCave.errMsg,
1474 <                  "SimSetup error: If you use a constant pressure\n"
1475 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1476 <        painCave.isFatal = 1;
1477 <        simError();
1478 <      }    
1479 <          
1480 <      if( globals->haveTauThermostat() )
1481 <        myNPTf->setTauThermostat( globals->getTauThermostat() );
1482 <      else{
1483 <        sprintf( painCave.errMsg,
1484 <         "SimSetup error: If you use an NPT\n"
1485 <                   "    ensemble, you must set tauThermostat.\n");
1486 <        painCave.isFatal = 1;
1487 <        simError();
1488 <      }
1489 <          
1490 <      if( globals->haveTauBarostat() )
1491 <        myNPTf->setTauBarostat( globals->getTauBarostat() );
1492 <      else{
1493 <        sprintf( painCave.errMsg,
1494 <                  "SimSetup error: If you use an NPT\n"
1495 <                  "    ensemble, you must set tauBarostat.\n");
1496 <        painCave.isFatal = 1;
1497 <        simError();
1498 <      }
1499 <      break;
1500 <      
1501 <    case NPTim_ENS:
1502 <      if (globals->haveZconstraints()){
1503 <        setupZConstraint(info[k]);
1504 <        myNPTim = new ZConstraint<NPTim<RealIntegrator> >( &(info[k]), the_ff );
1505 <      }
1506 <      else
1507 <        myNPTim = new NPTim<RealIntegrator>( &(info[k]), the_ff );
1391 >        myNPTi->setTargetTemp(globals->getTargetTemp());
1392  
1393 <        myNPTim->setTargetTemp( globals->getTargetTemp());
1394 <          
1395 <      if (globals->haveTargetPressure())
1396 <        myNPTim->setTargetPressure(globals->getTargetPressure());
1513 <      else {
1514 <        sprintf( painCave.errMsg,
1393 >        if (globals->haveTargetPressure())
1394 >          myNPTi->setTargetPressure(globals->getTargetPressure());
1395 >        else{
1396 >          sprintf(painCave.errMsg,
1397                    "SimSetup error: If you use a constant pressure\n"
1398                    "    ensemble, you must set targetPressure in the BASS file.\n");
1399 <        painCave.isFatal = 1;
1400 <        simError();
1401 <      }
1402 <          
1403 <      if( globals->haveTauThermostat() )
1404 <        myNPTim->setTauThermostat( globals->getTauThermostat() );
1405 <      else{
1406 <        sprintf( painCave.errMsg,
1399 >          painCave.isFatal = 1;
1400 >          simError();
1401 >        }
1402 >
1403 >        if (globals->haveTauThermostat())
1404 >          myNPTi->setTauThermostat(globals->getTauThermostat());
1405 >        else{
1406 >          sprintf(painCave.errMsg,
1407                    "SimSetup error: If you use an NPT\n"
1408                    "    ensemble, you must set tauThermostat.\n");
1409 <        painCave.isFatal = 1;
1410 <        simError();
1411 <      }
1530 <          
1531 <      if( globals->haveTauBarostat() )
1532 <        myNPTim->setTauBarostat( globals->getTauBarostat() );
1533 <      else{
1534 <        sprintf( painCave.errMsg,
1535 <                   "SimSetup error: If you use an NPT\n"
1536 <                   "    ensemble, you must set tauBarostat.\n");
1537 <        painCave.isFatal = 1;
1538 <        simError();
1539 <      }
1540 <      break;
1541 <      
1542 <    case NPTfm_ENS:
1543 <      if (globals->haveZconstraints()){
1544 <        setupZConstraint(info[k]);
1545 <        myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >( &(info[k]), the_ff );
1546 <      }
1547 <      else
1548 <        myNPTfm = new NPTfm<RealIntegrator>( &(info[k]), the_ff );
1409 >          painCave.isFatal = 1;
1410 >          simError();
1411 >        }
1412  
1413 <      myNPTfm->setTargetTemp( globals->getTargetTemp());
1413 >        if (globals->haveTauBarostat())
1414 >          myNPTi->setTauBarostat(globals->getTauBarostat());
1415 >        else{
1416 >          sprintf(painCave.errMsg,
1417 >                  "SimSetup error: If you use an NPT\n"
1418 >                  "    ensemble, you must set tauBarostat.\n");
1419 >          painCave.isFatal = 1;
1420 >          simError();
1421 >        }
1422 >        break;
1423  
1424 <      if (globals->haveTargetPressure())
1425 <        myNPTfm->setTargetPressure(globals->getTargetPressure());
1426 <      else {
1427 <        sprintf( painCave.errMsg,
1424 >      case NPTf_ENS:
1425 >        if (globals->haveZconstraints()){
1426 >          setupZConstraint(info[k]);
1427 >          myNPTf = new ZConstraint<NPTf<RealIntegrator> >(&(info[k]), the_ff);
1428 >        }
1429 >        else
1430 >          myNPTf = new NPTf<RealIntegrator>(&(info[k]), the_ff);
1431 >
1432 >        myNPTf->setTargetTemp(globals->getTargetTemp());
1433 >
1434 >        if (globals->haveTargetPressure())
1435 >          myNPTf->setTargetPressure(globals->getTargetPressure());
1436 >        else{
1437 >          sprintf(painCave.errMsg,
1438                    "SimSetup error: If you use a constant pressure\n"
1439                    "    ensemble, you must set targetPressure in the BASS file.\n");
1440 <        painCave.isFatal = 1;
1441 <        simError();
1442 <      }
1440 >          painCave.isFatal = 1;
1441 >          simError();
1442 >        }    
1443  
1444 <      if( globals->haveTauThermostat() )
1445 <        myNPTfm->setTauThermostat( globals->getTauThermostat() );
1446 <      else{
1447 <        sprintf( painCave.errMsg,
1444 >        if (globals->haveTauThermostat())
1445 >          myNPTf->setTauThermostat(globals->getTauThermostat());
1446 >        else{
1447 >          sprintf(painCave.errMsg,
1448                    "SimSetup error: If you use an NPT\n"
1449                    "    ensemble, you must set tauThermostat.\n");
1450 <        painCave.isFatal = 1;
1451 <        simError();
1452 <      }
1450 >          painCave.isFatal = 1;
1451 >          simError();
1452 >        }
1453  
1454 <      if( globals->haveTauBarostat() )
1455 <        myNPTfm->setTauBarostat( globals->getTauBarostat() );
1456 <      else{
1457 <        sprintf( painCave.errMsg,
1454 >        if (globals->haveTauBarostat())
1455 >          myNPTf->setTauBarostat(globals->getTauBarostat());
1456 >        else{
1457 >          sprintf(painCave.errMsg,
1458                    "SimSetup error: If you use an NPT\n"
1459                    "    ensemble, you must set tauBarostat.\n");
1460 +          painCave.isFatal = 1;
1461 +          simError();
1462 +        }
1463 +        break;
1464 +
1465 +      case NPTim_ENS:
1466 +        if (globals->haveZconstraints()){
1467 +          setupZConstraint(info[k]);
1468 +          myNPTim = new ZConstraint<NPTim<RealIntegrator> >(&(info[k]), the_ff);
1469 +        }
1470 +        else
1471 +          myNPTim = new NPTim<RealIntegrator>(&(info[k]), the_ff);
1472 +
1473 +        myNPTim->setTargetTemp(globals->getTargetTemp());
1474 +
1475 +        if (globals->haveTargetPressure())
1476 +          myNPTim->setTargetPressure(globals->getTargetPressure());
1477 +        else{
1478 +          sprintf(painCave.errMsg,
1479 +                  "SimSetup error: If you use a constant pressure\n"
1480 +                  "    ensemble, you must set targetPressure in the BASS file.\n");
1481 +          painCave.isFatal = 1;
1482 +          simError();
1483 +        }
1484 +
1485 +        if (globals->haveTauThermostat())
1486 +          myNPTim->setTauThermostat(globals->getTauThermostat());
1487 +        else{
1488 +          sprintf(painCave.errMsg,
1489 +                  "SimSetup error: If you use an NPT\n"
1490 +                  "    ensemble, you must set tauThermostat.\n");
1491 +          painCave.isFatal = 1;
1492 +          simError();
1493 +        }
1494 +
1495 +        if (globals->haveTauBarostat())
1496 +          myNPTim->setTauBarostat(globals->getTauBarostat());
1497 +        else{
1498 +          sprintf(painCave.errMsg,
1499 +                  "SimSetup error: If you use an NPT\n"
1500 +                  "    ensemble, you must set tauBarostat.\n");
1501 +          painCave.isFatal = 1;
1502 +          simError();
1503 +        }
1504 +        break;
1505 +
1506 +      case NPTfm_ENS:
1507 +        if (globals->haveZconstraints()){
1508 +          setupZConstraint(info[k]);
1509 +          myNPTfm = new ZConstraint<NPTfm<RealIntegrator> >(&(info[k]), the_ff);
1510 +        }
1511 +        else
1512 +          myNPTfm = new NPTfm<RealIntegrator>(&(info[k]), the_ff);
1513 +
1514 +        myNPTfm->setTargetTemp(globals->getTargetTemp());
1515 +
1516 +        if (globals->haveTargetPressure())
1517 +          myNPTfm->setTargetPressure(globals->getTargetPressure());
1518 +        else{
1519 +          sprintf(painCave.errMsg,
1520 +                  "SimSetup error: If you use a constant pressure\n"
1521 +                  "    ensemble, you must set targetPressure in the BASS file.\n");
1522 +          painCave.isFatal = 1;
1523 +          simError();
1524 +        }
1525 +
1526 +        if (globals->haveTauThermostat())
1527 +          myNPTfm->setTauThermostat(globals->getTauThermostat());
1528 +        else{
1529 +          sprintf(painCave.errMsg,
1530 +                  "SimSetup error: If you use an NPT\n"
1531 +                  "    ensemble, you must set tauThermostat.\n");
1532 +          painCave.isFatal = 1;
1533 +          simError();
1534 +        }
1535 +
1536 +        if (globals->haveTauBarostat())
1537 +          myNPTfm->setTauBarostat(globals->getTauBarostat());
1538 +        else{
1539 +          sprintf(painCave.errMsg,
1540 +                  "SimSetup error: If you use an NPT\n"
1541 +                  "    ensemble, you must set tauBarostat.\n");
1542 +          painCave.isFatal = 1;
1543 +          simError();
1544 +        }
1545 +        break;
1546 +
1547 +      default:
1548 +        sprintf(painCave.errMsg,
1549 +                "SimSetup Error. Unrecognized ensemble in case statement.\n");
1550          painCave.isFatal = 1;
1551          simError();
1580      }
1581      break;
1582      
1583    default:
1584      sprintf( painCave.errMsg,
1585                 "SimSetup Error. Unrecognized ensemble in case statement.\n");
1586      painCave.isFatal = 1;
1587      simError();
1552      }
1553    }
1554   }
1555  
1556 < void SimSetup::initFortran( void ){
1593 <
1556 > void SimSetup::initFortran(void){
1557    info[0].refreshSim();
1558 <  
1559 <  if( !strcmp( info[0].mixingRule, "standard") ){
1560 <    the_ff->initForceField( LB_MIXING_RULE );
1558 >
1559 >  if (!strcmp(info[0].mixingRule, "standard")){
1560 >    the_ff->initForceField(LB_MIXING_RULE);
1561    }
1562 <  else if( !strcmp( info[0].mixingRule, "explicit") ){
1563 <    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1562 >  else if (!strcmp(info[0].mixingRule, "explicit")){
1563 >    the_ff->initForceField(EXPLICIT_MIXING_RULE);
1564    }
1565    else{
1566 <    sprintf( painCave.errMsg,
1567 <       "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1605 <       info[0].mixingRule );
1566 >    sprintf(painCave.errMsg, "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1567 >            info[0].mixingRule);
1568      painCave.isFatal = 1;
1569      simError();
1570    }
1571  
1572  
1573   #ifdef IS_MPI
1574 <  strcpy( checkPointMsg,
1613 <    "Successfully intialized the mixingRule for Fortran." );
1574 >  strcpy(checkPointMsg, "Successfully intialized the mixingRule for Fortran.");
1575    MPIcheckPoint();
1576   #endif // is_mpi
1616
1577   }
1578  
1579 < void SimSetup::setupZConstraint(SimInfo& theInfo)
1620 < {
1579 > void SimSetup::setupZConstraint(SimInfo& theInfo){
1580    int nZConstraints;
1581    ZconStamp** zconStamp;
1582  
1583 <  if(globals->haveZconstraintTime()){  
1625 <    
1583 >  if (globals->haveZconstraintTime()){
1584      //add sample time of z-constraint  into SimInfo's property list                    
1585      DoubleData* zconsTimeProp = new DoubleData();
1586      zconsTimeProp->setID(ZCONSTIME_ID);
# Line 1630 | Line 1588 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1588      theInfo.addProperty(zconsTimeProp);
1589    }
1590    else{
1591 <    sprintf( painCave.errMsg,
1592 <       "ZConstraint error: If you use an ZConstraint\n"
1593 <       " , you must set sample time.\n");
1591 >    sprintf(painCave.errMsg,
1592 >            "ZConstraint error: If you use an ZConstraint\n"
1593 >            " , you must set sample time.\n");
1594      painCave.isFatal = 1;
1595 <    simError();      
1595 >    simError();
1596    }
1597  
1598    //push zconsTol into siminfo, if user does not specify
1599    //value for zconsTol, a default value will be used
1600    DoubleData* zconsTol = new DoubleData();
1601    zconsTol->setID(ZCONSTOL_ID);
1602 <  if(globals->haveZconsTol()){
1602 >  if (globals->haveZconsTol()){
1603      zconsTol->setData(globals->getZconsTol());
1604    }
1605    else{
1606 <  double defaultZConsTol = 0.01;
1607 <    sprintf( painCave.errMsg,
1608 <       "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1609 <       " , default value %f is used.\n", defaultZConsTol);
1606 >    double defaultZConsTol = 0.01;
1607 >    sprintf(painCave.errMsg,
1608 >            "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1609 >            " , default value %f is used.\n",
1610 >            defaultZConsTol);
1611      painCave.isFatal = 0;
1612      simError();      
1613  
# Line 1657 | Line 1616 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1616    theInfo.addProperty(zconsTol);
1617  
1618    //set Force Substraction Policy
1619 <  StringData* zconsForcePolicy =  new StringData();
1619 >  StringData* zconsForcePolicy = new StringData();
1620    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1621 <  
1622 <  if(globals->haveZconsForcePolicy()){
1621 >
1622 >  if (globals->haveZconsForcePolicy()){
1623      zconsForcePolicy->setData(globals->getZconsForcePolicy());
1624 <  }  
1624 >  }
1625    else{
1626 <     sprintf( painCave.errMsg,
1627 <             "ZConstraint Warning: User does not set force substraction policy, "
1628 <             "average force substraction policy is used\n");
1629 <     painCave.isFatal = 0;
1630 <     simError();
1631 <     zconsForcePolicy->setData("BYNUMBER");
1626 >    sprintf(painCave.errMsg,
1627 >            "ZConstraint Warning: User does not set force substraction policy, "
1628 >            "average force substraction policy is used\n");
1629 >    painCave.isFatal = 0;
1630 >    simError();
1631 >    zconsForcePolicy->setData("BYNUMBER");
1632    }
1633 <
1634 < theInfo.addProperty(zconsForcePolicy);
1635 <
1633 >
1634 >  theInfo.addProperty(zconsForcePolicy);
1635 >
1636    //Determine the name of ouput file and add it into SimInfo's property list
1637    //Be careful, do not use inFileName, since it is a pointer which
1638    //point to a string at master node, and slave nodes do not contain that string
1639 <  
1639 >
1640    string zconsOutput(theInfo.finalName);
1641 <  
1641 >
1642    zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1643 <  
1643 >
1644    StringData* zconsFilename = new StringData();
1645    zconsFilename->setID(ZCONSFILENAME_ID);
1646    zconsFilename->setData(zconsOutput);
1647 <  
1647 >
1648    theInfo.addProperty(zconsFilename);
1649 <  
1649 >
1650    //setup index, pos and other parameters of z-constraint molecules
1651    nZConstraints = globals->getNzConstraints();
1652    theInfo.nZconstraints = nZConstraints;
# Line 1697 | Line 1656 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1656  
1657    ZConsParaData* zconsParaData = new ZConsParaData();
1658    zconsParaData->setID(ZCONSPARADATA_ID);
1659 <
1660 <  for(int i = 0; i < nZConstraints; i++){
1659 >
1660 >  for (int i = 0; i < nZConstraints; i++){
1661      tempParaItem.havingZPos = zconStamp[i]->haveZpos();
1662      tempParaItem.zPos = zconStamp[i]->getZpos();
1663      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
# Line 1709 | Line 1668 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1668  
1669    //sort the parameters by index of molecules
1670    zconsParaData->sortByIndex();
1671 <  
1671 >
1672    //push data into siminfo, therefore, we can retrieve later
1673    theInfo.addProperty(zconsParaData);
1715      
1674   }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines