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 1187 by chrisfen, Sat May 22 18:16:18 2004 UTC

# Line 1 | Line 1
1   #include <algorithm>
2 < #include <cstdlib>
2 > #include <stdlib.h>
3   #include <iostream>
4 < #include <cmath>
4 > #include <math.h>
5   #include <string>
6 < #include <sys/time.h>
7 <
6 > #include <sprng.h>
7   #include "SimSetup.hpp"
8   #include "ReadWrite.hpp"
9   #include "parse_me.h"
10   #include "Integrator.hpp"
11   #include "simError.h"
12 + #include "RigidBody.hpp"
13 + //#include "ConjugateMinimizer.hpp"
14 + #include "OOPSEMinimizer.hpp"
15  
16   #ifdef IS_MPI
17   #include "mpiBASS.h"
# Line 22 | Line 24
24   #define NVT_ENS        1
25   #define NPTi_ENS       2
26   #define NPTf_ENS       3
27 < #define NPTim_ENS      4
26 < #define NPTfm_ENS      5
27 > #define NPTxyz_ENS     4
28  
28 – #define FF_DUFF 0
29 – #define FF_LJ   1
30 – #define FF_EAM  2
29  
30 + #define FF_DUFF  0
31 + #define FF_LJ    1
32 + #define FF_EAM   2
33 + #define FF_H2O   3
34 +
35   using namespace std;
36  
37 + /**
38 + * Check whether dividend is divisble by divisor or not
39 + */
40 + bool isDivisible(double dividend, double divisor){
41 +  double tolerance = 0.000001;
42 +  double quotient;
43 +  double diff;
44 +  int intQuotient;
45 +  
46 +  quotient = dividend / divisor;
47 +
48 +  if (quotient < 0)
49 +    quotient = -quotient;
50 +
51 +  intQuotient = int (quotient + tolerance);
52 +
53 +  diff = fabs(fabs(dividend) - intQuotient  * fabs(divisor));
54 +
55 +  if (diff <= tolerance)
56 +    return true;
57 +  else
58 +    return false;  
59 + }
60 +
61   SimSetup::SimSetup(){
62    
63 +  initSuspend = false;
64    isInfoArray = 0;
65    nInfo = 1;
66 <  
66 >
67    stamps = new MakeStamps();
68    globals = new Globals();
69 <  
70 <  
69 >
70 >
71   #ifdef IS_MPI
72 <  strcpy( checkPointMsg, "SimSetup creation successful" );
72 >  strcpy(checkPointMsg, "SimSetup creation successful");
73    MPIcheckPoint();
74   #endif // IS_MPI
75   }
# Line 51 | Line 79 | SimSetup::~SimSetup(){
79    delete globals;
80   }
81  
82 < void SimSetup::setSimInfo( SimInfo* the_info, int theNinfo ) {
83 <    info = the_info;
84 <    nInfo = theNinfo;
85 <    isInfoArray = 1;
82 > void SimSetup::setSimInfo(SimInfo* the_info, int theNinfo){
83 >  info = the_info;
84 >  nInfo = theNinfo;
85 >  isInfoArray = 1;
86 >  initSuspend = true;
87   }
88  
89  
90 < void SimSetup::parseFile( char* fileName ){
62 <
90 > void SimSetup::parseFile(char* fileName){
91   #ifdef IS_MPI
92 <  if( worldRank == 0 ){
92 >  if (worldRank == 0){
93   #endif // is_mpi
94 <    
94 >
95      inFileName = fileName;
96 <    set_interface_stamps( stamps, globals );
97 <    
96 >    set_interface_stamps(stamps, globals);
97 >
98   #ifdef IS_MPI
99      mpiEventInit();
100   #endif
101  
102 <    yacc_BASS( fileName );
102 >    yacc_BASS(fileName);
103  
104   #ifdef IS_MPI
105      throwMPIEvent(NULL);
106    }
107 <  else receiveParse();
107 >  else{
108 >    receiveParse();
109 >  }
110   #endif
111  
112   }
113  
114   #ifdef IS_MPI
115   void SimSetup::receiveParse(void){
116 <
117 <    set_interface_stamps( stamps, globals );
118 <    mpiEventInit();
119 <    MPIcheckPoint();
90 <    mpiEventLoop();
91 <
116 >  set_interface_stamps(stamps, globals);
117 >  mpiEventInit();
118 >  MPIcheckPoint();
119 >  mpiEventLoop();
120   }
121  
122   #endif // is_mpi
123  
124   void SimSetup::createSim(void){
125  
98 –  int i, j, k, globalAtomIndex;
99 –  
126    // gather all of the information from the Bass file
127  
128    gatherInfo();
# Line 111 | Line 137 | void SimSetup::createSim(void){
137  
138    // initialize the system coordinates
139  
140 <  if( !isInfoArray ) initSystemCoords();  
140 >  if ( !initSuspend ){
141 >    initSystemCoords();
142  
143 +    if( !(globals->getUseInitTime()) )
144 +      info[0].currentTime = 0.0;
145 +  }  
146 +
147    // make the output filenames
148  
149    makeOutNames();
119 –  
120 –  // make the integrator
150    
122 –  makeIntegrator();
123 –  
151   #ifdef IS_MPI
152    mpiSim->mpiRefresh();
153   #endif
# Line 129 | Line 156 | void SimSetup::createSim(void){
156  
157    initFortran();
158  
159 +  if (globals->haveMinimizer())
160 +    // make minimizer
161 +    makeMinimizer();
162 +  else
163 +    // make the integrator
164 +    makeIntegrator();
165  
133 –
166   }
167  
168  
169 < void SimSetup::makeMolecules( void ){
170 <
171 <  int k,l;
172 <  int i, j, exI, exJ, tempEx, stampID, atomOffset, excludeOffset;
169 > void SimSetup::makeMolecules(void){
170 >  int i, j, k;
171 >  int exI, exJ, exK, exL, slI, slJ;
172 >  int tempI, tempJ, tempK, tempL;
173 >  int molI;
174 >  int stampID, atomOffset, rbOffset;
175    molInit molInfo;
176    DirectionalAtom* dAtom;
177 +  RigidBody* myRB;
178 +  StuntDouble* mySD;
179    LinkedAssign* extras;
180    LinkedAssign* current_extra;
181    AtomStamp* currentAtom;
182    BondStamp* currentBond;
183    BendStamp* currentBend;
184    TorsionStamp* currentTorsion;
185 +  RigidBodyStamp* currentRigidBody;
186 +  CutoffGroupStamp* currentCutoffGroup;
187 +  CutoffGroup* myCutoffGroup;
188 +  int nCutoffGroups;// number of cutoff group of a molecule defined in mdl file
189 +  set<int> cutoffAtomSet; //atoms belong to  cutoffgroup defined at mdl file
190  
191    bond_pair* theBonds;
192    bend_set* theBends;
193    torsion_set* theTorsions;
194  
195 <  
195 >  set<int> skipList;
196 >
197 >  double phi, theta, psi;
198 >  char* molName;
199 >  char rbName[100];
200 >
201    //init the forceField paramters
202  
203    the_ff->readParams();
204  
159 –  
205    // init the atoms
206  
207 <  double ux, uy, uz, u, uSqr;
208 <  
209 <  for(k=0; k<nInfo; k++){
210 <    
166 <    the_ff->setSimInfo( &(info[k]) );
207 >  int nMembers, nNew, rb1, rb2;
208 >
209 >  for (k = 0; k < nInfo; k++){
210 >    the_ff->setSimInfo(&(info[k]));
211  
212      atomOffset = 0;
213 <    excludeOffset = 0;
214 <    for(i=0; i<info[k].n_mol; i++){
171 <    
213 >
214 >    for (i = 0; i < info[k].n_mol; i++){
215        stampID = info[k].molecules[i].getStampID();
216 +      molName = comp_stamps[stampID]->getID();
217  
218 <      molInfo.nAtoms    = comp_stamps[stampID]->getNAtoms();
219 <      molInfo.nBonds    = comp_stamps[stampID]->getNBonds();
220 <      molInfo.nBends    = comp_stamps[stampID]->getNBends();
218 >      molInfo.nAtoms = comp_stamps[stampID]->getNAtoms();
219 >      molInfo.nBonds = comp_stamps[stampID]->getNBonds();
220 >      molInfo.nBends = comp_stamps[stampID]->getNBends();
221        molInfo.nTorsions = comp_stamps[stampID]->getNTorsions();
222 <      molInfo.nExcludes = molInfo.nBonds + molInfo.nBends + molInfo.nTorsions;
222 >      molInfo.nRigidBodies = comp_stamps[stampID]->getNRigidBodies();
223 >
224 >      nCutoffGroups = comp_stamps[stampID]->getNCutoffGroups();
225        
226        molInfo.myAtoms = &(info[k].atoms[atomOffset]);
181 –      molInfo.myExcludes = &(info[k].excludes[excludeOffset]);
182 –      molInfo.myBonds = new Bond*[molInfo.nBonds];
183 –      molInfo.myBends = new Bend*[molInfo.nBends];
184 –      molInfo.myTorsions = new Torsion*[molInfo.nTorsions];
227  
228 +      if (molInfo.nBonds > 0)
229 +        molInfo.myBonds = new Bond*[molInfo.nBonds];
230 +      else
231 +        molInfo.myBonds = NULL;
232 +
233 +      if (molInfo.nBends > 0)
234 +        molInfo.myBends = new Bend*[molInfo.nBends];
235 +      else
236 +        molInfo.myBends = NULL;
237 +
238 +      if (molInfo.nTorsions > 0)
239 +        molInfo.myTorsions = new Torsion *[molInfo.nTorsions];
240 +      else
241 +        molInfo.myTorsions = NULL;
242 +
243        theBonds = new bond_pair[molInfo.nBonds];
244        theBends = new bend_set[molInfo.nBends];
245        theTorsions = new torsion_set[molInfo.nTorsions];
246 <    
246 >      
247        // make the Atoms
248 <    
249 <      for(j=0; j<molInfo.nAtoms; j++){
250 <  
251 <  currentAtom = comp_stamps[stampID]->getAtom( j );
252 <  if( currentAtom->haveOrientation() ){
253 <    
254 <    dAtom = new DirectionalAtom( (j + atomOffset),
255 <               info[k].getConfiguration() );
256 <    info[k].n_oriented++;
257 <    molInfo.myAtoms[j] = dAtom;
258 <    
259 <    ux = currentAtom->getOrntX();
260 <    uy = currentAtom->getOrntY();
261 <    uz = currentAtom->getOrntZ();
262 <    
263 <    uSqr = (ux * ux) + (uy * uy) + (uz * uz);
264 <    
265 <    u = sqrt( uSqr );
266 <    ux = ux / u;
267 <    uy = uy / u;
268 <    uz = uz / u;
269 <    
270 <    dAtom->setSUx( ux );
271 <    dAtom->setSUy( uy );
272 <    dAtom->setSUz( uz );
273 <  }
274 <  else{
275 <    molInfo.myAtoms[j] = new GeneralAtom( (j + atomOffset),
219 <            info[k].getConfiguration() );
220 <  }
221 <  molInfo.myAtoms[j]->setType( currentAtom->getType() );
222 <    
248 >
249 >      for (j = 0; j < molInfo.nAtoms; j++){
250 >        currentAtom = comp_stamps[stampID]->getAtom(j);
251 >
252 >        if (currentAtom->haveOrientation()){
253 >          dAtom = new DirectionalAtom((j + atomOffset),
254 >                                      info[k].getConfiguration());
255 >          info[k].n_oriented++;
256 >          molInfo.myAtoms[j] = dAtom;
257 >
258 >          // Directional Atoms have standard unit vectors which are oriented
259 >          // in space using the three Euler angles.  We assume the standard
260 >          // unit vector was originally along the z axis below.
261 >
262 >          phi = currentAtom->getEulerPhi() * M_PI / 180.0;
263 >          theta = currentAtom->getEulerTheta() * M_PI / 180.0;
264 >          psi = currentAtom->getEulerPsi()* M_PI / 180.0;
265 >
266 >          dAtom->setUnitFrameFromEuler(phi, theta, psi);
267 >            
268 >        }
269 >        else{
270 >
271 >          molInfo.myAtoms[j] = new Atom((j + atomOffset), info[k].getConfiguration());
272 >
273 >        }
274 >
275 >        molInfo.myAtoms[j]->setType(currentAtom->getType());
276   #ifdef IS_MPI
277 <      
278 <  molInfo.myAtoms[j]->setGlobalIndex( globalIndex[j+atomOffset] );
279 <      
277 >
278 >        molInfo.myAtoms[j]->setGlobalIndex(globalAtomIndex[j + atomOffset]);
279 >
280   #endif // is_mpi
281        }
282 <    
283 <    // make the bonds
284 <      for(j=0; j<molInfo.nBonds; j++){
285 <      
286 <  currentBond = comp_stamps[stampID]->getBond( j );
287 <  theBonds[j].a = currentBond->getA() + atomOffset;
288 <  theBonds[j].b = currentBond->getB() + atomOffset;
289 <  
290 <  exI = theBonds[j].a;
291 <  exJ = theBonds[j].b;
292 <  
293 <  // exclude_I must always be the smaller of the pair
294 <  if( exI > exJ ){
295 <    tempEx = exI;
296 <    exI = exJ;
297 <    exJ = tempEx;
298 <  }
299 < #ifdef IS_MPI
300 <  tempEx = exI;
301 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
302 <  tempEx = exJ;
250 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
251 <  
252 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
253 < #else  // isn't MPI
254 <  
255 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
256 < #endif  //is_mpi
257 <      }
258 <      excludeOffset += molInfo.nBonds;
259 <      
282 >
283 >      // make the bonds
284 >      for (j = 0; j < molInfo.nBonds; j++){
285 >        currentBond = comp_stamps[stampID]->getBond(j);
286 >        theBonds[j].a = currentBond->getA() + atomOffset;
287 >        theBonds[j].b = currentBond->getB() + atomOffset;
288 >
289 >        tempI = theBonds[j].a;
290 >        tempJ = theBonds[j].b;
291 >
292 > #ifdef IS_MPI
293 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
294 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
295 > #else
296 >        exI = tempI + 1;
297 >        exJ = tempJ + 1;
298 > #endif
299 >
300 >        info[k].excludes->addPair(exI, exJ);
301 >      }
302 >
303        //make the bends
304 <      for(j=0; j<molInfo.nBends; j++){
305 <  
306 <  currentBend = comp_stamps[stampID]->getBend( j );
307 <  theBends[j].a = currentBend->getA() + atomOffset;
308 <  theBends[j].b = currentBend->getB() + atomOffset;
309 <  theBends[j].c = currentBend->getC() + atomOffset;
310 <  
311 <  if( currentBend->haveExtras() ){
312 <    
313 <    extras = currentBend->getExtras();
314 <    current_extra = extras;
315 <    
316 <    while( current_extra != NULL ){
317 <      if( !strcmp( current_extra->getlhs(), "ghostVectorSource" )){
304 >      for (j = 0; j < molInfo.nBends; j++){
305 >        currentBend = comp_stamps[stampID]->getBend(j);
306 >        theBends[j].a = currentBend->getA() + atomOffset;
307 >        theBends[j].b = currentBend->getB() + atomOffset;
308 >        theBends[j].c = currentBend->getC() + atomOffset;
309 >
310 >        if (currentBend->haveExtras()){
311 >          extras = currentBend->getExtras();
312 >          current_extra = extras;
313 >
314 >          while (current_extra != NULL){
315 >            if (!strcmp(current_extra->getlhs(), "ghostVectorSource")){
316 >              switch (current_extra->getType()){
317 >                case 0:
318 >                  theBends[j].ghost = current_extra->getInt() + atomOffset;
319 >                  theBends[j].isGhost = 1;
320 >                  break;
321 >
322 >                case 1:
323 >                  theBends[j].ghost = (int) current_extra->getDouble() +
324 >                                      atomOffset;
325 >                  theBends[j].isGhost = 1;
326 >                  break;
327 >
328 >                default:
329 >                  sprintf(painCave.errMsg,
330 >                          "SimSetup Error: ghostVectorSource was neither a "
331 >                          "double nor an int.\n"
332 >                          "-->Bend[%d] in %s\n",
333 >                          j, comp_stamps[stampID]->getID());
334 >                  painCave.isFatal = 1;
335 >                  simError();
336 >              }
337 >            }
338 >            else{
339 >              sprintf(painCave.errMsg,
340 >                      "SimSetup Error: unhandled bend assignment:\n"
341 >                      "    -->%s in Bend[%d] in %s\n",
342 >                      current_extra->getlhs(), j, comp_stamps[stampID]->getID());
343 >              painCave.isFatal = 1;
344 >              simError();
345 >            }
346 >
347 >            current_extra = current_extra->getNext();
348 >          }
349 >        }
350 >
351 >        if (theBends[j].isGhost) {
352 >          
353 >          tempI = theBends[j].a;
354 >          tempJ = theBends[j].b;
355 >          
356 > #ifdef IS_MPI
357 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
358 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
359 > #else
360 >          exI = tempI + 1;
361 >          exJ = tempJ + 1;
362 > #endif          
363 >          info[k].excludes->addPair(exI, exJ);
364 >
365 >        } else {
366 >
367 >          tempI = theBends[j].a;
368 >          tempJ = theBends[j].b;
369 >          tempK = theBends[j].c;
370 >          
371 > #ifdef IS_MPI
372 >          exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
373 >          exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
374 >          exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
375 > #else
376 >          exI = tempI + 1;
377 >          exJ = tempJ + 1;
378 >          exK = tempK + 1;
379 > #endif
380 >          
381 >          info[k].excludes->addPair(exI, exK);
382 >          info[k].excludes->addPair(exI, exJ);
383 >          info[k].excludes->addPair(exJ, exK);
384 >        }
385 >      }
386 >
387 >      for (j = 0; j < molInfo.nTorsions; j++){
388 >        currentTorsion = comp_stamps[stampID]->getTorsion(j);
389 >        theTorsions[j].a = currentTorsion->getA() + atomOffset;
390 >        theTorsions[j].b = currentTorsion->getB() + atomOffset;
391 >        theTorsions[j].c = currentTorsion->getC() + atomOffset;
392 >        theTorsions[j].d = currentTorsion->getD() + atomOffset;
393 >
394 >        tempI = theTorsions[j].a;      
395 >        tempJ = theTorsions[j].b;
396 >        tempK = theTorsions[j].c;
397 >        tempL = theTorsions[j].d;
398 >
399 > #ifdef IS_MPI
400 >        exI = info[k].atoms[tempI]->getGlobalIndex() + 1;
401 >        exJ = info[k].atoms[tempJ]->getGlobalIndex() + 1;
402 >        exK = info[k].atoms[tempK]->getGlobalIndex() + 1;
403 >        exL = info[k].atoms[tempL]->getGlobalIndex() + 1;
404 > #else
405 >        exI = tempI + 1;
406 >        exJ = tempJ + 1;
407 >        exK = tempK + 1;
408 >        exL = tempL + 1;
409 > #endif
410 >
411 >        info[k].excludes->addPair(exI, exJ);
412 >        info[k].excludes->addPair(exI, exK);
413 >        info[k].excludes->addPair(exI, exL);        
414 >        info[k].excludes->addPair(exJ, exK);
415 >        info[k].excludes->addPair(exJ, exL);
416 >        info[k].excludes->addPair(exK, exL);
417 >      }
418 >
419 >      
420 >      molInfo.myRigidBodies.clear();
421 >      
422 >      for (j = 0; j < molInfo.nRigidBodies; j++){
423 >
424 >        currentRigidBody = comp_stamps[stampID]->getRigidBody(j);
425 >        nMembers = currentRigidBody->getNMembers();
426 >
427 >        // Create the Rigid Body:
428 >
429 >        myRB = new RigidBody();
430 >
431 >        sprintf(rbName,"%s_RB_%d", molName, j);
432 >        myRB->setType(rbName);
433          
434 <        switch( current_extra->getType() ){
435 <    
436 <        case 0:
437 <    theBends[j].ghost =
438 <      current_extra->getInt() + atomOffset;
439 <    theBends[j].isGhost = 1;
440 <    break;
441 <    
442 <        case 1:
443 <    theBends[j].ghost =
444 <      (int)current_extra->getDouble() + atomOffset;
445 <    theBends[j].isGhost = 1;
446 <    break;
447 <    
448 <        default:
449 <    sprintf( painCave.errMsg,
450 <       "SimSetup Error: ghostVectorSource was neither a "
451 <       "double nor an int.\n"
452 <       "-->Bend[%d] in %s\n",
453 <       j, comp_stamps[stampID]->getID() );
454 <    painCave.isFatal = 1;
455 <    simError();
434 >        for (rb1 = 0; rb1 < nMembers; rb1++) {
435 >
436 >          // molI is atom numbering inside this molecule
437 >          molI = currentRigidBody->getMember(rb1);    
438 >
439 >          // tempI is atom numbering on local processor
440 >          tempI = molI + atomOffset;
441 >
442 >          // currentAtom is the AtomStamp (which we need for
443 >          // rigid body reference positions)
444 >          currentAtom = comp_stamps[stampID]->getAtom(molI);
445 >
446 >          // When we add to the rigid body, add the atom itself and
447 >          // the stamp info:
448 >
449 >          myRB->addAtom(info[k].atoms[tempI], currentAtom);
450 >          
451 >          // Add this atom to the Skip List for the integrators
452 > #ifdef IS_MPI
453 >          slI = info[k].atoms[tempI]->getGlobalIndex();
454 > #else
455 >          slI = tempI;
456 > #endif
457 >          skipList.insert(slI);
458 >          
459          }
460 +        
461 +        for(rb1 = 0; rb1 < nMembers - 1; rb1++) {
462 +          for(rb2 = rb1+1; rb2 < nMembers; rb2++) {
463 +            
464 +            tempI = currentRigidBody->getMember(rb1);
465 +            tempJ = currentRigidBody->getMember(rb2);
466 +            
467 +            // Some explanation is required here.
468 +            // Fortran indexing starts at 1, while c indexing starts at 0
469 +            // Also, in parallel computations, the GlobalIndex is
470 +            // used for the exclude list:
471 +            
472 + #ifdef IS_MPI
473 +            exI = molInfo.myAtoms[tempI]->getGlobalIndex() + 1;
474 +            exJ = molInfo.myAtoms[tempJ]->getGlobalIndex() + 1;
475 + #else
476 +            exI = molInfo.myAtoms[tempI]->getIndex() + 1;
477 +            exJ = molInfo.myAtoms[tempJ]->getIndex() + 1;
478 + #endif
479 +            
480 +            info[k].excludes->addPair(exI, exJ);
481 +            
482 +          }
483 +        }
484 +
485 +        molInfo.myRigidBodies.push_back(myRB);
486 +        info[k].rigidBodies.push_back(myRB);
487        }
488        
489 <      else{
489 >
490 >      //create cutoff group for molecule
491 >
492 >      cutoffAtomSet.clear();
493 >      molInfo.myCutoffGroups.clear();
494 >      
495 >      for (j = 0; j < nCutoffGroups; j++){
496 >
497 >        currentCutoffGroup = comp_stamps[stampID]->getCutoffGroup(j);
498 >        nMembers = currentCutoffGroup->getNMembers();
499 >
500 >        myCutoffGroup = new CutoffGroup();
501          
502 <        sprintf( painCave.errMsg,
503 <           "SimSetup Error: unhandled bend assignment:\n"
504 <           "    -->%s in Bend[%d] in %s\n",
505 <           current_extra->getlhs(),
506 <           j, comp_stamps[stampID]->getID() );
507 <        painCave.isFatal = 1;
508 <        simError();
502 >        for (int cg = 0; cg < nMembers; cg++) {
503 >
504 >          // molI is atom numbering inside this molecule
505 >          molI = currentCutoffGroup->getMember(cg);    
506 >
507 >          // tempI is atom numbering on local processor
508 >          tempI = molI + atomOffset;
509 >          
510 >          myCutoffGroup->addAtom(info[k].atoms[tempI]);          
511 >
512 >          cutoffAtomSet.insert(tempI);
513 >        }
514 >
515 >        molInfo.myCutoffGroups.push_back(myCutoffGroup);
516 >      }//end for (j = 0; j < molInfo.nCutoffGroups; j++)
517 >
518 >      //creat a cutoff group for every atom  in current molecule which does not belong to cutoffgroup defined at mdl file
519 >
520 >      for(j = 0; j < molInfo.nAtoms; j++){
521 >
522 >        if(cutoffAtomSet.find(molInfo.myAtoms[j]->getIndex()) == cutoffAtomSet.end()){
523 >          myCutoffGroup = new CutoffGroup();
524 >          myCutoffGroup->addAtom(molInfo.myAtoms[j]);
525 >          molInfo.myCutoffGroups.push_back(myCutoffGroup);
526 >        }
527 >          
528        }
529 +
530 +              
531 +
532 +
533 +      // After this is all set up, scan through the atoms to
534 +      // see if they can be added to the integrableObjects:
535 +
536 +      molInfo.myIntegrableObjects.clear();
537        
538 <      current_extra = current_extra->getNext();
539 <    }
540 <  }
315 <  
316 <  if( !theBends[j].isGhost ){
317 <    
318 <    exI = theBends[j].a;
319 <    exJ = theBends[j].c;
320 <  }
321 <  else{
322 <    
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 <  }
538 >
539 >      for (j = 0; j < molInfo.nAtoms; j++){
540 >
541   #ifdef IS_MPI
542 <  tempEx = exI;
543 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
544 <  tempEx = exJ;
545 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
546 <      
547 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
548 < #else  // isn't MPI
549 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
550 < #endif  //is_mpi
542 >        slJ = molInfo.myAtoms[j]->getGlobalIndex();
543 > #else
544 >        slJ = j+atomOffset;
545 > #endif
546 >
547 >        // if they aren't on the skip list, then they can be integrated
548 >
549 >        if (skipList.find(slJ) == skipList.end()) {
550 >          mySD = (StuntDouble *) molInfo.myAtoms[j];
551 >          info[k].integrableObjects.push_back(mySD);
552 >          molInfo.myIntegrableObjects.push_back(mySD);
553 >        }
554        }
555 <      excludeOffset += molInfo.nBends;
556 <      
557 <      for(j=0; j<molInfo.nTorsions; j++){
558 <  
559 <  currentTorsion = comp_stamps[stampID]->getTorsion( j );
560 <  theTorsions[j].a = currentTorsion->getA() + atomOffset;
561 <  theTorsions[j].b = currentTorsion->getB() + atomOffset;
351 <  theTorsions[j].c = currentTorsion->getC() + atomOffset;
352 <  theTorsions[j].d = currentTorsion->getD() + atomOffset;
353 <  
354 <  exI = theTorsions[j].a;
355 <  exJ = theTorsions[j].d;
356 <  
357 <  // exclude_I must always be the smaller of the pair
358 <  if( exI > exJ ){
359 <    tempEx = exI;
360 <    exI = exJ;
361 <    exJ = tempEx;
362 <  }
363 < #ifdef IS_MPI
364 <  tempEx = exI;
365 <  exI = info[k].atoms[tempEx]->getGlobalIndex() + 1;
366 <  tempEx = exJ;
367 <  exJ = info[k].atoms[tempEx]->getGlobalIndex() + 1;
368 <  
369 <  info[k].excludes[j+excludeOffset]->setPair( exI, exJ );
370 < #else  // isn't MPI
371 <  info[k].excludes[j+excludeOffset]->setPair( (exI+1), (exJ+1) );
372 < #endif  //is_mpi
555 >
556 >      // all rigid bodies are integrated:
557 >
558 >      for (j = 0; j < molInfo.nRigidBodies; j++) {
559 >        mySD = (StuntDouble *) molInfo.myRigidBodies[j];
560 >        info[k].integrableObjects.push_back(mySD);      
561 >        molInfo.myIntegrableObjects.push_back(mySD);
562        }
563 <      excludeOffset += molInfo.nTorsions;
563 >    
564        
376 –      
565        // send the arrays off to the forceField for init.
566        
567 <      the_ff->initializeAtoms( molInfo.nAtoms, molInfo.myAtoms );
568 <      the_ff->initializeBonds( molInfo.nBonds, molInfo.myBonds, theBonds );
569 <      the_ff->initializeBends( molInfo.nBends, molInfo.myBends, theBends );
570 <      the_ff->initializeTorsions( molInfo.nTorsions, molInfo.myTorsions, theTorsions );
571 <      
384 <      
385 <      info[k].molecules[i].initialize( molInfo );
567 >      the_ff->initializeAtoms(molInfo.nAtoms, molInfo.myAtoms);
568 >      the_ff->initializeBonds(molInfo.nBonds, molInfo.myBonds, theBonds);
569 >      the_ff->initializeBends(molInfo.nBends, molInfo.myBends, theBends);
570 >      the_ff->initializeTorsions(molInfo.nTorsions, molInfo.myTorsions,
571 >                                 theTorsions);
572  
573 <      
573 >      info[k].molecules[i].initialize(molInfo);
574 >
575 >
576        atomOffset += molInfo.nAtoms;
577        delete[] theBonds;
578        delete[] theBends;
579        delete[] theTorsions;
580 <    }
580 >    }    
581    }
582 <  
582 >
583   #ifdef IS_MPI
584 <  sprintf( checkPointMsg, "all molecules initialized succesfully" );
584 >  sprintf(checkPointMsg, "all molecules initialized succesfully");
585    MPIcheckPoint();
586   #endif // is_mpi
399 –  
400 –  // clean up the forcefield
587  
402 –  the_ff->calcRcut();
403 –  the_ff->cleanMe();
404 –  
588   }
589  
590 < void SimSetup::initFromBass( void ){
408 <
590 > void SimSetup::initFromBass(void){
591    int i, j, k;
592    int n_cells;
593    double cellx, celly, cellz;
# Line 419 | Line 601 | void SimSetup::initFromBass( void ){
601    vel[1] = 0.0;
602    vel[2] = 0.0;
603  
604 <  temp1 = (double)tot_nmol / 4.0;
605 <  temp2 = pow( temp1, ( 1.0 / 3.0 ) );
606 <  temp3 = ceil( temp2 );
604 >  temp1 = (double) tot_nmol / 4.0;
605 >  temp2 = pow(temp1, (1.0 / 3.0));
606 >  temp3 = ceil(temp2);
607  
608 <  have_extra =0;
609 <  if( temp2 < temp3 ){ // we have a non-complete lattice
610 <    have_extra =1;
608 >  have_extra = 0;
609 >  if (temp2 < temp3){
610 >    // we have a non-complete lattice
611 >    have_extra = 1;
612  
613 <    n_cells = (int)temp3 - 1;
613 >    n_cells = (int) temp3 - 1;
614      cellx = info[0].boxL[0] / temp3;
615      celly = info[0].boxL[1] / temp3;
616      cellz = info[0].boxL[2] / temp3;
617 <    n_extra = tot_nmol - ( 4 * n_cells * n_cells * n_cells );
618 <    temp1 = ((double)n_extra) / ( pow( temp3, 3.0 ) - pow( n_cells, 3.0 ) );
619 <    n_per_extra = (int)ceil( temp1 );
617 >    n_extra = tot_nmol - (4 * n_cells * n_cells * n_cells);
618 >    temp1 = ((double) n_extra) / (pow(temp3, 3.0) - pow(n_cells, 3.0));
619 >    n_per_extra = (int) ceil(temp1);
620  
621 <    if( n_per_extra > 4){
622 <      sprintf( painCave.errMsg,
623 <         "SimSetup error. There has been an error in constructing"
624 <         " the non-complete lattice.\n" );
621 >    if (n_per_extra > 4){
622 >      sprintf(painCave.errMsg,
623 >              "SimSetup error. There has been an error in constructing"
624 >              " the non-complete lattice.\n");
625        painCave.isFatal = 1;
626        simError();
627      }
628    }
629    else{
630 <    n_cells = (int)temp3;
630 >    n_cells = (int) temp3;
631      cellx = info[0].boxL[0] / temp3;
632      celly = info[0].boxL[1] / temp3;
633      cellz = info[0].boxL[2] / temp3;
# Line 455 | Line 638 | void SimSetup::initFromBass( void ){
638    current_comp = 0;
639    current_atom_ndx = 0;
640  
641 <  for( i=0; i < n_cells ; i++ ){
642 <    for( j=0; j < n_cells; j++ ){
643 <      for( k=0; k < n_cells; k++ ){
641 >  for (i = 0; i < n_cells ; i++){
642 >    for (j = 0; j < n_cells; j++){
643 >      for (k = 0; k < n_cells; k++){
644 >        makeElement(i * cellx, j * celly, k * cellz);
645  
646 <  makeElement( i * cellx,
463 <         j * celly,
464 <         k * cellz );
646 >        makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly, k * cellz);
647  
648 <  makeElement( i * cellx + 0.5 * cellx,
467 <         j * celly + 0.5 * celly,
468 <         k * cellz );
648 >        makeElement(i * cellx, j * celly + 0.5 * celly, k * cellz + 0.5 * cellz);
649  
650 <  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 );
650 >        makeElement(i * cellx + 0.5 * cellx, j * celly, k * cellz + 0.5 * cellz);
651        }
652      }
653    }
654  
655 <  if( have_extra ){
655 >  if (have_extra){
656      done = 0;
657  
658      int start_ndx;
659 <    for( i=0; i < (n_cells+1) && !done; i++ ){
660 <      for( j=0; j < (n_cells+1) && !done; j++ ){
659 >    for (i = 0; i < (n_cells + 1) && !done; i++){
660 >      for (j = 0; j < (n_cells + 1) && !done; j++){
661 >        if (i < n_cells){
662 >          if (j < n_cells){
663 >            start_ndx = n_cells;
664 >          }
665 >          else
666 >            start_ndx = 0;
667 >        }
668 >        else
669 >          start_ndx = 0;
670  
671 <  if( i < n_cells ){
671 >        for (k = start_ndx; k < (n_cells + 1) && !done; k++){
672 >          makeElement(i * cellx, j * celly, k * cellz);
673 >          done = (current_mol >= tot_nmol);
674  
675 <    if( j < n_cells ){
676 <      start_ndx = n_cells;
677 <    }
678 <    else start_ndx = 0;
679 <  }
495 <  else start_ndx = 0;
675 >          if (!done && n_per_extra > 1){
676 >            makeElement(i * cellx + 0.5 * cellx, j * celly + 0.5 * celly,
677 >                        k * cellz);
678 >            done = (current_mol >= tot_nmol);
679 >          }
680  
681 <  for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
681 >          if (!done && n_per_extra > 2){
682 >            makeElement(i * cellx, j * celly + 0.5 * celly,
683 >                        k * cellz + 0.5 * cellz);
684 >            done = (current_mol >= tot_nmol);
685 >          }
686  
687 <    makeElement( i * cellx,
688 <           j * celly,
689 <           k * cellz );
690 <    done = ( current_mol >= tot_nmol );
691 <
692 <    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 <    }
517 <
518 <    if( !done && n_per_extra > 3){
519 <      makeElement( i * cellx + 0.5 * cellx,
520 <       j * celly,
521 <       k * cellz + 0.5 * cellz );
522 <      done = ( current_mol >= tot_nmol );
523 <    }
524 <  }
687 >          if (!done && n_per_extra > 3){
688 >            makeElement(i * cellx + 0.5 * cellx, j * celly,
689 >                        k * cellz + 0.5 * cellz);
690 >            done = (current_mol >= tot_nmol);
691 >          }
692 >        }
693        }
694      }
695    }
696  
697 <  for( i=0; i<info[0].n_atoms; i++ ){
698 <    info[0].atoms[i]->setVel( vel );
697 >  for (i = 0; i < info[0].n_atoms; i++){
698 >    info[0].atoms[i]->setVel(vel);
699    }
700   }
701  
702 < void SimSetup::makeElement( double x, double y, double z ){
535 <
702 > void SimSetup::makeElement(double x, double y, double z){
703    int k;
704    AtomStamp* current_atom;
705    DirectionalAtom* dAtom;
706    double rotMat[3][3];
707    double pos[3];
708  
709 <  for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
710 <
711 <    current_atom = comp_stamps[current_comp]->getAtom( k );
712 <    if( !current_atom->havePosition() ){
713 <      sprintf( painCave.errMsg,
714 <         "SimSetup:initFromBass error.\n"
715 <         "\tComponent %s, atom %s does not have a position specified.\n"
716 <         "\tThe initialization routine is unable to give a start"
717 <         " position.\n",
551 <         comp_stamps[current_comp]->getID(),
552 <         current_atom->getType() );
709 >  for (k = 0; k < comp_stamps[current_comp]->getNAtoms(); k++){
710 >    current_atom = comp_stamps[current_comp]->getAtom(k);
711 >    if (!current_atom->havePosition()){
712 >      sprintf(painCave.errMsg,
713 >              "SimSetup:initFromBass error.\n"
714 >              "\tComponent %s, atom %s does not have a position specified.\n"
715 >              "\tThe initialization routine is unable to give a start"
716 >              " position.\n",
717 >              comp_stamps[current_comp]->getID(), current_atom->getType());
718        painCave.isFatal = 1;
719        simError();
720      }
721 <    
721 >
722      pos[0] = x + current_atom->getPosX();
723      pos[1] = y + current_atom->getPosY();
724      pos[2] = z + current_atom->getPosZ();
560 –    
561 –    info[0].atoms[current_atom_ndx]->setPos( pos );
725  
726 <    if( info[0].atoms[current_atom_ndx]->isDirectional() ){
726 >    info[0].atoms[current_atom_ndx]->setPos(pos);
727  
728 <      dAtom = (DirectionalAtom *)info[0].atoms[current_atom_ndx];
728 >    if (info[0].atoms[current_atom_ndx]->isDirectional()){
729 >      dAtom = (DirectionalAtom *) info[0].atoms[current_atom_ndx];
730  
731        rotMat[0][0] = 1.0;
732        rotMat[0][1] = 0.0;
# Line 576 | Line 740 | void SimSetup::makeElement( double x, double y, double
740        rotMat[2][1] = 0.0;
741        rotMat[2][2] = 1.0;
742  
743 <      dAtom->setA( rotMat );
743 >      dAtom->setA(rotMat);
744      }
745  
746      current_atom_ndx++;
# Line 585 | Line 749 | void SimSetup::makeElement( double x, double y, double
749    current_mol++;
750    current_comp_mol++;
751  
752 <  if( current_comp_mol >= components_nmol[current_comp] ){
589 <
752 >  if (current_comp_mol >= components_nmol[current_comp]){
753      current_comp_mol = 0;
754      current_comp++;
755    }
756   }
757  
758  
759 < void SimSetup::gatherInfo( void ){
760 <  int i,j,k;
759 > void SimSetup::gatherInfo(void){
760 >  int i;
761  
762    ensembleCase = -1;
763    ffCase = -1;
764  
765    // set the easy ones first
766  
767 <  for( i=0; i<nInfo; i++){
767 >  for (i = 0; i < nInfo; i++){
768      info[i].target_temp = globals->getTargetTemp();
769      info[i].dt = globals->getDt();
770      info[i].run_time = globals->getRunTime();
# Line 611 | Line 774 | void SimSetup::gatherInfo( void ){
774  
775    // get the forceField
776  
777 <  strcpy( force_field, globals->getForceField() );
777 >  strcpy(force_field, globals->getForceField());
778  
779 <  if( !strcasecmp( force_field, "DUFF" )) ffCase = FF_DUFF;
780 <  else if( !strcasecmp( force_field, "LJ" )) ffCase = FF_LJ;
781 <  else if( !strcasecmp( force_field, "EAM" )) ffCase = FF_EAM;
779 >  if (!strcasecmp(force_field, "DUFF")){
780 >    ffCase = FF_DUFF;
781 >  }
782 >  else if (!strcasecmp(force_field, "LJ")){
783 >    ffCase = FF_LJ;
784 >  }
785 >  else if (!strcasecmp(force_field, "EAM")){
786 >    ffCase = FF_EAM;
787 >  }
788 >  else if (!strcasecmp(force_field, "WATER")){
789 >    ffCase = FF_H2O;
790 >  }
791    else{
792 <    sprintf( painCave.errMsg,
793 <       "SimSetup Error. Unrecognized force field -> %s\n",
794 <       force_field );
795 <    painCave.isFatal = 1;
624 <    simError();
792 >    sprintf(painCave.errMsg, "SimSetup Error. Unrecognized force field -> %s\n",
793 >            force_field);
794 >         painCave.isFatal = 1;
795 >         simError();
796    }
797  
798 <  // get the ensemble
798 >    // get the ensemble
799  
800 <  strcpy( ensemble, globals->getEnsemble() );
800 >  strcpy(ensemble, globals->getEnsemble());
801  
802 <  if( !strcasecmp( ensemble, "NVE" ))      ensembleCase = NVE_ENS;
803 <  else if( !strcasecmp( ensemble, "NVT" )) ensembleCase = NVT_ENS;
804 <  else if( !strcasecmp( ensemble, "NPTi" ) || !strcasecmp( ensemble, "NPT") )
802 >  if (!strcasecmp(ensemble, "NVE")){
803 >    ensembleCase = NVE_ENS;
804 >  }
805 >  else if (!strcasecmp(ensemble, "NVT")){
806 >    ensembleCase = NVT_ENS;
807 >  }
808 >  else if (!strcasecmp(ensemble, "NPTi") || !strcasecmp(ensemble, "NPT")){
809      ensembleCase = NPTi_ENS;
810 <  else if( !strcasecmp( ensemble, "NPTf" )) ensembleCase = NPTf_ENS;
811 <  else if( !strcasecmp( ensemble, "NPTim" )) ensembleCase = NPTim_ENS;
812 <  else if( !strcasecmp( ensemble, "NPTfm" )) ensembleCase = NPTfm_ENS;
810 >  }
811 >  else if (!strcasecmp(ensemble, "NPTf")){
812 >    ensembleCase = NPTf_ENS;
813 >  }
814 >  else if (!strcasecmp(ensemble, "NPTxyz")){
815 >    ensembleCase = NPTxyz_ENS;
816 >  }
817    else{
818 <    sprintf( painCave.errMsg,
819 <       "SimSetup Warning. Unrecognized Ensemble -> %s, "
820 <             "reverting to NVE for this simulation.\n",
821 <       ensemble );
822 <    painCave.isFatal = 0;
823 <    simError();
824 <    strcpy( ensemble, "NVE" );
825 <    ensembleCase = NVE_ENS;
818 >    sprintf(painCave.errMsg,
819 >            "SimSetup Warning. Unrecognized Ensemble -> %s \n"
820 >            "\treverting to NVE for this simulation.\n",
821 >            ensemble);
822 >         painCave.isFatal = 0;
823 >         simError();
824 >         strcpy(ensemble, "NVE");
825 >         ensembleCase = NVE_ENS;
826    }  
648 –  
649 –  for(i=0; i<nInfo; i++){
650 –    
651 –    strcpy( info[i].ensemble, ensemble );
827  
828 +  for (i = 0; i < nInfo; i++){
829 +    strcpy(info[i].ensemble, ensemble);
830 +
831      // get the mixing rule
832  
833 <    strcpy( info[i].mixingRule, globals->getMixingRule() );
833 >    strcpy(info[i].mixingRule, globals->getMixingRule());
834      info[i].usePBC = globals->getPBC();
835    }
836 <  
836 >
837    // get the components and calculate the tot_nMol and indvidual n_mol
838 <
838 >
839    the_components = globals->getComponents();
840    components_nmol = new int[n_components];
841  
842  
843 <  if( !globals->haveNMol() ){
843 >  if (!globals->haveNMol()){
844      // we don't have the total number of molecules, so we assume it is
845      // given in each component
846  
847      tot_nmol = 0;
848 <    for( i=0; i<n_components; i++ ){
849 <
850 <      if( !the_components[i]->haveNMol() ){
851 <  // we have a problem
852 <  sprintf( painCave.errMsg,
853 <     "SimSetup Error. No global NMol or component NMol"
854 <     " given. Cannot calculate the number of atoms.\n" );
855 <  painCave.isFatal = 1;
678 <  simError();
848 >    for (i = 0; i < n_components; i++){
849 >      if (!the_components[i]->haveNMol()){
850 >        // we have a problem
851 >        sprintf(painCave.errMsg,
852 >                "SimSetup Error. No global NMol or component NMol given.\n"
853 >                "\tCannot calculate the number of atoms.\n");
854 >        painCave.isFatal = 1;
855 >        simError();
856        }
857  
858        tot_nmol += the_components[i]->getNMol();
# Line 683 | Line 860 | void SimSetup::gatherInfo( void ){
860      }
861    }
862    else{
863 <    sprintf( painCave.errMsg,
864 <       "SimSetup error.\n"
865 <       "\tSorry, the ability to specify total"
866 <       " nMols and then give molfractions in the components\n"
867 <       "\tis not currently supported."
868 <       " Please give nMol in the components.\n" );
863 >    sprintf(painCave.errMsg,
864 >            "SimSetup error.\n"
865 >            "\tSorry, the ability to specify total"
866 >            " nMols and then give molfractions in the components\n"
867 >            "\tis not currently supported."
868 >            " Please give nMol in the components.\n");
869      painCave.isFatal = 1;
870      simError();
871    }
872  
873 +  //check whether sample time, status time, thermal time and reset time are divisble by dt
874 +  if (globals->haveSampleTime() && !isDivisible(globals->getSampleTime(), globals->getDt())){
875 +    sprintf(painCave.errMsg,
876 +            "Sample time is not divisible by dt.\n"
877 +            "\tThis will result in samples that are not uniformly\n"
878 +            "\tdistributed in time.  If this is a problem, change\n"
879 +            "\tyour sampleTime variable.\n");
880 +    painCave.isFatal = 0;
881 +    simError();    
882 +  }
883 +
884 +  if (globals->haveStatusTime() && !isDivisible(globals->getStatusTime(), globals->getDt())){
885 +    sprintf(painCave.errMsg,
886 +            "Status time is not divisible by dt.\n"
887 +            "\tThis will result in status reports that are not uniformly\n"
888 +            "\tdistributed in time.  If this is a problem, change \n"
889 +            "\tyour statusTime variable.\n");
890 +    painCave.isFatal = 0;
891 +    simError();    
892 +  }
893 +
894 +  if (globals->haveThermalTime() && !isDivisible(globals->getThermalTime(), globals->getDt())){
895 +    sprintf(painCave.errMsg,
896 +            "Thermal time is not divisible by dt.\n"
897 +            "\tThis will result in thermalizations that are not uniformly\n"
898 +            "\tdistributed in time.  If this is a problem, change \n"
899 +            "\tyour thermalTime variable.\n");
900 +    painCave.isFatal = 0;
901 +    simError();    
902 +  }  
903 +
904 +  if (globals->haveResetTime() && !isDivisible(globals->getResetTime(), globals->getDt())){
905 +    sprintf(painCave.errMsg,
906 +            "Reset time is not divisible by dt.\n"
907 +            "\tThis will result in integrator resets that are not uniformly\n"
908 +            "\tdistributed in time.  If this is a problem, change\n"
909 +            "\tyour resetTime variable.\n");
910 +    painCave.isFatal = 0;
911 +    simError();    
912 +  }
913 +
914    // set the status, sample, and thermal kick times
697 –  
698 –  for(i=0; i<nInfo; i++){
915  
916 <    if( globals->haveSampleTime() ){
916 >  for (i = 0; i < nInfo; i++){
917 >    if (globals->haveSampleTime()){
918        info[i].sampleTime = globals->getSampleTime();
919        info[i].statusTime = info[i].sampleTime;
703 –      info[i].thermalTime = info[i].sampleTime;
920      }
921      else{
922        info[i].sampleTime = globals->getRunTime();
923        info[i].statusTime = info[i].sampleTime;
708 –      info[i].thermalTime = info[i].sampleTime;
924      }
925 <    
926 <    if( globals->haveStatusTime() ){
925 >
926 >    if (globals->haveStatusTime()){
927        info[i].statusTime = globals->getStatusTime();
928      }
929 <    
930 <    if( globals->haveThermalTime() ){
929 >
930 >    if (globals->haveThermalTime()){
931        info[i].thermalTime = globals->getThermalTime();
932 +    } else {
933 +      info[i].thermalTime = globals->getRunTime();
934      }
935  
936 <    // check for the temperature set flag
936 >    info[i].resetIntegrator = 0;
937 >    if( globals->haveResetTime() ){
938 >      info[i].resetTime = globals->getResetTime();
939 >      info[i].resetIntegrator = 1;
940 >    }
941  
942 <    if( globals->haveTempSet() ) info[i].setTemp = globals->getTempSet();
942 >    // check for the temperature set flag
943      
944 <    // get some of the tricky things that may still be in the globals
945 <    
946 <    double boxVector[3];
947 <    if( globals->haveBox() ){
948 <      boxVector[0] = globals->getBox();
949 <      boxVector[1] = globals->getBox();
950 <      boxVector[2] = globals->getBox();
951 <      
952 <      info[i].setBox( boxVector );
953 <    }
954 <    else if( globals->haveDensity() ){
955 <      
956 <      double vol;
957 <      vol = (double)tot_nmol / globals->getDensity();
958 <      boxVector[0] = pow( vol, ( 1.0 / 3.0 ) );
959 <      boxVector[1] = boxVector[0];
960 <      boxVector[2] = boxVector[0];
740 <      
741 <      info[i].setBox( boxVector );
742 <  }
743 <    else{
744 <      if( !globals->haveBoxX() ){
745 <  sprintf( painCave.errMsg,
746 <     "SimSetup error, no periodic BoxX size given.\n" );
747 <  painCave.isFatal = 1;
748 <  simError();
944 >    if (globals->haveTempSet())
945 >      info[i].setTemp = globals->getTempSet();
946 >
947 >    // check for the extended State init
948 >
949 >    info[i].useInitXSstate = globals->getUseInitXSstate();
950 >    info[i].orthoTolerance = globals->getOrthoBoxTolerance();
951 >
952 >    // check for thermodynamic integration
953 >    if (globals->getUseThermInt()) {
954 >      if (globals->haveThermIntLambda() && globals->haveThermIntK()) {
955 >        info[i].useThermInt = globals->getUseThermInt();
956 >        info[i].thermIntLambda = globals->getThermIntLambda();
957 >        info[i].thermIntK = globals->getThermIntK();
958 >        
959 >        Restraints *myRestraint = new Restraints(tot_nmol, info[i].thermIntLambda, info[i].thermIntK);
960 >        info[i].restraint = myRestraint;
961        }
962 <      boxVector[0] = globals->getBoxX();
963 <      
964 <      if( !globals->haveBoxY() ){
965 <  sprintf( painCave.errMsg,
966 <     "SimSetup error, no periodic BoxY size given.\n" );
967 <  painCave.isFatal = 1;
968 <  simError();
962 >      else {
963 >        sprintf(painCave.errMsg,
964 >                "SimSetup Error:\n"
965 >                "\tKeyword useThermInt was set to 'true' but\n"
966 >                "\tthermodynamicIntegrationLambda (and/or\n"
967 >                "\tthermodynamicIntegrationK) was not specified.\n"
968 >                "\tPlease provide a lambda value and k value in your .bass file.\n");
969 >        painCave.isFatal = 1;
970 >        simError();    
971        }
758 –      boxVector[1] = globals->getBoxY();
759 –      
760 –      if( !globals->haveBoxZ() ){
761 –  sprintf( painCave.errMsg,
762 –     "SimSetup error, no periodic BoxZ size given.\n" );
763 –  painCave.isFatal = 1;
764 –  simError();
765 –      }
766 –      boxVector[2] = globals->getBoxZ();
767 –      
768 –      info[i].setBox( boxVector );
972      }
973 +    else if(globals->haveThermIntLambda() || globals->haveThermIntK()){
974 +        sprintf(painCave.errMsg,
975 +                "SimSetup Warning: If you want to use Thermodynamic\n"
976 +                "\tIntegration, set useThermInt to 'true' in your .bass file.\n"
977 +                "\tThe useThermInt keyword is 'false' by default, so your\n"
978 +                "\tlambda and/or k values are being ignored.\n");
979 +        painCave.isFatal = 0;
980 +        simError();  
981 +    }
982    }
983 <
983 >  
984 >  //setup seed for random number generator
985    int seedValue;
773 –  struct timeval now_time_val;
774 –  struct timezone time_zone;
986  
987 <  if(globals->haveSeed()){
987 >  if (globals->haveSeed()){
988      seedValue = globals->getSeed();
989 <  }
989 >
990 >    if(seedValue / 1E9 == 0){
991 >      sprintf(painCave.errMsg,
992 >              "Seed for sprng library should contain at least 9 digits\n"
993 >              "OOPSE will generate a seed for user\n");
994 >      painCave.isFatal = 0;
995 >      simError();
996 >
997 >      //using seed generated by system instead of invalid seed set by user
998 > #ifndef IS_MPI
999 >      seedValue = make_sprng_seed();
1000 > #else
1001 >      if (worldRank == 0){
1002 >        seedValue = make_sprng_seed();
1003 >      }
1004 >      MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
1005 > #endif      
1006 >    }
1007 >  }//end of if branch of globals->haveSeed()
1008    else{
1009 +    
1010   #ifndef IS_MPI
1011 <    gettimeofday(&now_time_val, &time_zone);  // get the time now  
782 <    seedValue = (int) now_time_val.tv_sec;      //  convert to epoch time
1011 >    seedValue = make_sprng_seed();
1012   #else
1013 <    if(worldRank == 0){
1014 <      gettimeofday(&now_time_val, &time_zone);  // get the time now
786 <      seedValue = (int) now_time_val.tv_sec;      //  convert to epoch time      
1013 >    if (worldRank == 0){
1014 >      seedValue = make_sprng_seed();
1015      }
1016 <     MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
1016 >    MPI_Bcast(&seedValue, 1, MPI_INT, 0, MPI_COMM_WORLD);  
1017   #endif
1018 <  }
1018 >  }//end of globals->haveSeed()
1019  
1020 <  for(int i = 0; i < nInfo; i++){
1020 >  for (int i = 0; i < nInfo; i++){
1021      info[i].setSeed(seedValue);
1022    }
1023 <
1023 >  
1024   #ifdef IS_MPI
1025 <  strcpy( checkPointMsg, "Succesfully gathered all information from Bass\n" );
1025 >  strcpy(checkPointMsg, "Successfully gathered all information from Bass\n");
1026    MPIcheckPoint();
1027   #endif // is_mpi
800 –
1028   }
1029  
1030  
1031 < void SimSetup::finalInfoCheck( void ){
1031 > void SimSetup::finalInfoCheck(void){
1032    int index;
1033    int usesDipoles;
1034 +  int usesCharges;
1035    int i;
1036  
1037 <  for(i=0; i<nInfo; i++){
1037 >  for (i = 0; i < nInfo; i++){
1038      // check electrostatic parameters
1039 <    
1039 >
1040      index = 0;
1041      usesDipoles = 0;
1042 <    while( (index < info[i].n_atoms) && !usesDipoles ){
1042 >    while ((index < info[i].n_atoms) && !usesDipoles){
1043        usesDipoles = (info[i].atoms[index])->hasDipole();
1044        index++;
1045      }
1046 <    
1046 >    index = 0;
1047 >    usesCharges = 0;
1048 >    while ((index < info[i].n_atoms) && !usesCharges){
1049 >      usesCharges= (info[i].atoms[index])->hasCharge();
1050 >      index++;
1051 >    }
1052   #ifdef IS_MPI
1053      int myUse = usesDipoles;
1054 <    MPI_Allreduce( &myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD );
1054 >    MPI_Allreduce(&myUse, &usesDipoles, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
1055   #endif //is_mpi
1056 <    
1057 <    double theEcr, theEst;
1058 <  
1059 <    if (globals->getUseRF() ) {
1060 <      info[i].useReactionField = 1;
1056 >
1057 >    double theRcut, theRsw;
1058 >
1059 >    if (globals->haveRcut()) {
1060 >      theRcut = globals->getRcut();
1061 >
1062 >      if (globals->haveRsw())
1063 >        theRsw = globals->getRsw();
1064 >      else
1065 >        theRsw = theRcut;
1066        
1067 <      if( !globals->haveECR() ){
1068 <  sprintf( painCave.errMsg,
1069 <     "SimSetup Warning: using default value of 1/2 the smallest "
832 <     "box length for the electrostaticCutoffRadius.\n"
833 <     "I hope you have a very fast processor!\n");
834 <  painCave.isFatal = 0;
835 <  simError();
836 <  double smallest;
837 <  smallest = info[i].boxL[0];
838 <  if (info[i].boxL[1] <= smallest) smallest = info[i].boxL[1];
839 <  if (info[i].boxL[2] <= smallest) smallest = info[i].boxL[2];
840 <  theEcr = 0.5 * smallest;
841 <      } else {
842 <  theEcr = globals->getECR();
843 <      }
1067 >      info[i].setDefaultRcut(theRcut, theRsw);
1068 >
1069 >    } else {
1070        
1071 <      if( !globals->haveEST() ){
1072 <  sprintf( painCave.errMsg,
1073 <     "SimSetup Warning: using default value of 0.05 * the "
1074 <     "electrostaticCutoffRadius for the electrostaticSkinThickness\n"
1075 <     );
1076 <  painCave.isFatal = 0;
1077 <  simError();
852 <  theEst = 0.05 * theEcr;
853 <      } else {
854 <  theEst= globals->getEST();
855 <      }
1071 >      the_ff->calcRcut();
1072 >      theRcut = info[i].getRcut();
1073 >
1074 >      if (globals->haveRsw())
1075 >        theRsw = globals->getRsw();
1076 >      else
1077 >        theRsw = theRcut;
1078        
1079 <      info[i].setEcr( theEcr, theEst );
1079 >      info[i].setDefaultRcut(theRcut, theRsw);
1080 >    }
1081 >
1082 >    if (globals->getUseRF()){
1083 >      info[i].useReactionField = 1;
1084        
1085 <      if(!globals->haveDielectric() ){
1086 <  sprintf( painCave.errMsg,
1087 <     "SimSetup Error: You are trying to use Reaction Field without"
1088 <     "setting a dielectric constant!\n"
1089 <     );
1090 <  painCave.isFatal = 1;
1091 <  simError();
1085 >      if (!globals->haveRcut()){
1086 >        sprintf(painCave.errMsg,
1087 >                "SimSetup Warning: No value was set for the cutoffRadius.\n"
1088 >                "\tOOPSE will use a default value of 15.0 angstroms"
1089 >                "\tfor the cutoffRadius.\n");
1090 >        painCave.isFatal = 0;
1091 >        simError();
1092 >        theRcut = 15.0;
1093        }
1094 <      info[i].dielectric = globals->getDielectric();  
1095 <    }
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 );
1094 >      else{
1095 >        theRcut = globals->getRcut();
1096        }
903 –    }  
904 –  }
1097  
1098 +      if (!globals->haveRsw()){
1099 +        sprintf(painCave.errMsg,
1100 +                "SimSetup Warning: No value was set for switchingRadius.\n"
1101 +                "\tOOPSE will use a default value of\n"
1102 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
1103 +        painCave.isFatal = 0;
1104 +        simError();
1105 +        theRsw = 0.95 * theRcut;
1106 +      }
1107 +      else{
1108 +        theRsw = globals->getRsw();
1109 +      }
1110 +
1111 +      info[i].setDefaultRcut(theRcut, theRsw);
1112 +
1113 +      if (!globals->haveDielectric()){
1114 +        sprintf(painCave.errMsg,
1115 +                "SimSetup Error: No Dielectric constant was set.\n"
1116 +                "\tYou are trying to use Reaction Field without"
1117 +                "\tsetting a dielectric constant!\n");
1118 +        painCave.isFatal = 1;
1119 +        simError();
1120 +      }
1121 +      info[i].dielectric = globals->getDielectric();
1122 +    }
1123 +    else{
1124 +      if (usesDipoles || usesCharges){
1125 +
1126 +        if (!globals->haveRcut()){
1127 +          sprintf(painCave.errMsg,
1128 +                  "SimSetup Warning: No value was set for the cutoffRadius.\n"
1129 +                  "\tOOPSE will use a default value of 15.0 angstroms"
1130 +                  "\tfor the cutoffRadius.\n");
1131 +          painCave.isFatal = 0;
1132 +          simError();
1133 +          theRcut = 15.0;
1134 +      }
1135 +        else{
1136 +          theRcut = globals->getRcut();
1137 +        }
1138 +        
1139 +        if (!globals->haveRsw()){
1140 +          sprintf(painCave.errMsg,
1141 +                  "SimSetup Warning: No value was set for switchingRadius.\n"
1142 +                  "\tOOPSE will use a default value of\n"
1143 +                  "\t0.95 * cutoffRadius for the switchingRadius\n");
1144 +          painCave.isFatal = 0;
1145 +          simError();
1146 +          theRsw = 0.95 * theRcut;
1147 +        }
1148 +        else{
1149 +          theRsw = globals->getRsw();
1150 +        }
1151 +        
1152 +        info[i].setDefaultRcut(theRcut, theRsw);
1153 +        
1154 +      }
1155 +    }
1156 +  }
1157   #ifdef IS_MPI
1158 <  strcpy( checkPointMsg, "post processing checks out" );
1158 >  strcpy(checkPointMsg, "post processing checks out");
1159    MPIcheckPoint();
1160   #endif // is_mpi
1161  
1162 +  // clean up the forcefield
1163 +  the_ff->cleanMe();
1164   }
912 –
913 – void SimSetup::initSystemCoords( void ){
914 –  int i;
1165    
1166 + void SimSetup::initSystemCoords(void){
1167 +  int i;
1168 +
1169    char* inName;
1170  
1171 <  (info[0].getConfiguration())->createArrays( info[0].n_atoms );
1171 >  (info[0].getConfiguration())->createArrays(info[0].n_atoms);
1172  
1173 <  for(i=0; i<info[0].n_atoms; i++) info[0].atoms[i]->setCoords();
1174 <  
1175 <  if( globals->haveInitialConfig() ){
1176 <    
1173 >  for (i = 0; i < info[0].n_atoms; i++)
1174 >    info[0].atoms[i]->setCoords();
1175 >
1176 >  if (globals->haveInitialConfig()){
1177      InitializeFromFile* fileInit;
1178   #ifdef IS_MPI // is_mpi
1179 <    if( worldRank == 0 ){
1179 >    if (worldRank == 0){
1180   #endif //is_mpi
1181        inName = globals->getInitialConfig();
1182 <      double* tempDouble = new double[1000000];
930 <      fileInit = new InitializeFromFile( inName );
1182 >      fileInit = new InitializeFromFile(inName);
1183   #ifdef IS_MPI
1184 <    }else fileInit = new InitializeFromFile( NULL );
1184 >    }
1185 >    else
1186 >      fileInit = new InitializeFromFile(NULL);
1187   #endif
1188 <    fileInit->readInit( info ); // default velocities on
1189 <    
1188 >    fileInit->readInit(info); // default velocities on
1189 >
1190      delete fileInit;
1191    }
1192    else{
1193      
940 – #ifdef IS_MPI
941 –    
1194      // no init from bass
1195      
1196 <    sprintf( painCave.errMsg,
1197 <       "Cannot intialize a parallel simulation without an initial configuration file.\n" );
1198 <    painCave.isFatal;
1196 >    sprintf(painCave.errMsg,
1197 >            "Cannot intialize a simulation without an initial configuration file.\n");
1198 >    painCave.isFatal = 1;;
1199      simError();
1200      
949 – #else
950 –    
951 –    initFromBass();
952 –    
953 –    
954 – #endif
1201    }
1202 <  
1202 >
1203   #ifdef IS_MPI
1204 <  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
1204 >  strcpy(checkPointMsg, "Successfully read in the initial configuration");
1205    MPIcheckPoint();
1206   #endif // is_mpi
961 –  
1207   }
1208  
1209  
1210 < void SimSetup::makeOutNames( void ){
966 <  
1210 > void SimSetup::makeOutNames(void){
1211    int k;
1212  
969 –  
970 –  for(k=0; k<nInfo; k++){
1213  
1214 +  for (k = 0; k < nInfo; k++){
1215   #ifdef IS_MPI
1216 <    if( worldRank == 0 ){
1216 >    if (worldRank == 0){
1217   #endif // is_mpi
1218 <      
1219 <      if( globals->haveFinalConfig() ){
1220 <  strcpy( info[k].finalName, globals->getFinalConfig() );
1218 >
1219 >      if (globals->haveFinalConfig()){
1220 >        strcpy(info[k].finalName, globals->getFinalConfig());
1221        }
1222        else{
1223 <  strcpy( info[k].finalName, inFileName );
1224 <  char* endTest;
1225 <  int nameLength = strlen( info[k].finalName );
1226 <  endTest = &(info[k].finalName[nameLength - 5]);
1227 <  if( !strcmp( endTest, ".bass" ) ){
1228 <    strcpy( endTest, ".eor" );
1229 <  }
1230 <  else if( !strcmp( endTest, ".BASS" ) ){
1231 <    strcpy( endTest, ".eor" );
1232 <  }
1233 <  else{
1234 <    endTest = &(info[k].finalName[nameLength - 4]);
1235 <    if( !strcmp( endTest, ".bss" ) ){
1236 <      strcpy( endTest, ".eor" );
1237 <    }
1238 <    else if( !strcmp( endTest, ".mdl" ) ){
1239 <      strcpy( endTest, ".eor" );
1240 <    }
1241 <    else{
1242 <      strcat( info[k].finalName, ".eor" );
1243 <    }
1244 <  }
1223 >        strcpy(info[k].finalName, inFileName);
1224 >        char* endTest;
1225 >        int nameLength = strlen(info[k].finalName);
1226 >        endTest = &(info[k].finalName[nameLength - 5]);
1227 >        if (!strcmp(endTest, ".bass")){
1228 >          strcpy(endTest, ".eor");
1229 >        }
1230 >        else if (!strcmp(endTest, ".BASS")){
1231 >          strcpy(endTest, ".eor");
1232 >        }
1233 >        else{
1234 >          endTest = &(info[k].finalName[nameLength - 4]);
1235 >          if (!strcmp(endTest, ".bss")){
1236 >            strcpy(endTest, ".eor");
1237 >          }
1238 >          else if (!strcmp(endTest, ".mdl")){
1239 >            strcpy(endTest, ".eor");
1240 >          }
1241 >          else{
1242 >            strcat(info[k].finalName, ".eor");
1243 >          }
1244 >        }
1245        }
1246 <      
1246 >
1247        // make the sample and status out names
1248 <      
1249 <      strcpy( info[k].sampleName, inFileName );
1248 >
1249 >      strcpy(info[k].sampleName, inFileName);
1250        char* endTest;
1251 <      int nameLength = strlen( info[k].sampleName );
1251 >      int nameLength = strlen(info[k].sampleName);
1252        endTest = &(info[k].sampleName[nameLength - 5]);
1253 <      if( !strcmp( endTest, ".bass" ) ){
1254 <  strcpy( endTest, ".dump" );
1253 >      if (!strcmp(endTest, ".bass")){
1254 >        strcpy(endTest, ".dump");
1255        }
1256 <      else if( !strcmp( endTest, ".BASS" ) ){
1257 <  strcpy( endTest, ".dump" );
1256 >      else if (!strcmp(endTest, ".BASS")){
1257 >        strcpy(endTest, ".dump");
1258        }
1259        else{
1260 <  endTest = &(info[k].sampleName[nameLength - 4]);
1261 <  if( !strcmp( endTest, ".bss" ) ){
1262 <    strcpy( endTest, ".dump" );
1263 <  }
1264 <  else if( !strcmp( endTest, ".mdl" ) ){
1265 <    strcpy( endTest, ".dump" );
1266 <  }
1267 <  else{
1268 <    strcat( info[k].sampleName, ".dump" );
1269 <  }
1260 >        endTest = &(info[k].sampleName[nameLength - 4]);
1261 >        if (!strcmp(endTest, ".bss")){
1262 >          strcpy(endTest, ".dump");
1263 >        }
1264 >        else if (!strcmp(endTest, ".mdl")){
1265 >          strcpy(endTest, ".dump");
1266 >        }
1267 >        else{
1268 >          strcat(info[k].sampleName, ".dump");
1269 >        }
1270        }
1271 <      
1272 <      strcpy( info[k].statusName, inFileName );
1273 <      nameLength = strlen( info[k].statusName );
1271 >
1272 >      strcpy(info[k].statusName, inFileName);
1273 >      nameLength = strlen(info[k].statusName);
1274        endTest = &(info[k].statusName[nameLength - 5]);
1275 <      if( !strcmp( endTest, ".bass" ) ){
1276 <  strcpy( endTest, ".stat" );
1275 >      if (!strcmp(endTest, ".bass")){
1276 >        strcpy(endTest, ".stat");
1277        }
1278 <      else if( !strcmp( endTest, ".BASS" ) ){
1279 <  strcpy( endTest, ".stat" );
1278 >      else if (!strcmp(endTest, ".BASS")){
1279 >        strcpy(endTest, ".stat");
1280        }
1281        else{
1282 <  endTest = &(info[k].statusName[nameLength - 4]);
1283 <  if( !strcmp( endTest, ".bss" ) ){
1284 <    strcpy( endTest, ".stat" );
1285 <  }
1286 <  else if( !strcmp( endTest, ".mdl" ) ){
1287 <    strcpy( endTest, ".stat" );
1288 <  }
1289 <  else{
1290 <    strcat( info[k].statusName, ".stat" );
1291 <  }
1282 >        endTest = &(info[k].statusName[nameLength - 4]);
1283 >        if (!strcmp(endTest, ".bss")){
1284 >          strcpy(endTest, ".stat");
1285 >        }
1286 >        else if (!strcmp(endTest, ".mdl")){
1287 >          strcpy(endTest, ".stat");
1288 >        }
1289 >        else{
1290 >          strcat(info[k].statusName, ".stat");
1291 >        }
1292        }
1293 <      
1293 >
1294 >      strcpy(info[k].rawPotName, inFileName);
1295 >      nameLength = strlen(info[k].rawPotName);
1296 >      endTest = &(info[k].rawPotName[nameLength - 5]);
1297 >      if (!strcmp(endTest, ".bass")){
1298 >        strcpy(endTest, ".raw");
1299 >      }
1300 >      else if (!strcmp(endTest, ".BASS")){
1301 >        strcpy(endTest, ".raw");
1302 >      }
1303 >      else{
1304 >        endTest = &(info[k].rawPotName[nameLength - 4]);
1305 >        if (!strcmp(endTest, ".bss")){
1306 >          strcpy(endTest, ".raw");
1307 >        }
1308 >        else if (!strcmp(endTest, ".mdl")){
1309 >          strcpy(endTest, ".raw");
1310 >        }
1311 >        else{
1312 >          strcat(info[k].rawPotName, ".raw");
1313 >        }
1314 >      }
1315 >
1316   #ifdef IS_MPI
1317 +
1318      }
1319   #endif // is_mpi
1320    }
1321   }
1322  
1323  
1324 < void SimSetup::sysObjectsCreation( void ){
1325 <  
1326 <  int i,k;
1061 <  
1324 > void SimSetup::sysObjectsCreation(void){
1325 >  int i, k;
1326 >
1327    // create the forceField
1328  
1329    createFF();
# Line 1073 | Line 1338 | void SimSetup::sysObjectsCreation( void ){
1338  
1339   #ifdef IS_MPI
1340    // divide the molecules among the processors
1341 <  
1341 >
1342    mpiMolDivide();
1343   #endif //is_mpi
1344 <  
1344 >
1345    // create the atom and SRI arrays. Also initialize Molecule Stamp ID's
1346 <
1346 >
1347    makeSysArrays();
1348  
1349    // make and initialize the molecules (all but atomic coordinates)
1350 <
1350 >
1351    makeMolecules();
1352 <  
1353 <  for(k=0; k<nInfo; k++){
1352 >
1353 >  for (k = 0; k < nInfo; k++){
1354      info[k].identArray = new int[info[k].n_atoms];
1355 <    for(i=0; i<info[k].n_atoms; i++){
1355 >    for (i = 0; i < info[k].n_atoms; i++){
1356        info[k].identArray[i] = info[k].atoms[i]->getIdent();
1357      }
1358    }
1359   }
1360  
1361  
1362 < void SimSetup::createFF( void ){
1362 > void SimSetup::createFF(void){
1363 >  switch (ffCase){
1364 >    case FF_DUFF:
1365 >      the_ff = new DUFF();
1366 >      break;
1367  
1368 <  switch( ffCase ){
1368 >    case FF_LJ:
1369 >      the_ff = new LJFF();
1370 >      break;
1371  
1372 <  case FF_DUFF:
1373 <    the_ff = new DUFF();
1374 <    break;
1372 >    case FF_EAM:
1373 >      the_ff = new EAM_FF();
1374 >      break;
1375  
1376 <  case FF_LJ:
1377 <    the_ff = new LJFF();
1378 <    break;
1376 >    case FF_H2O:
1377 >      the_ff = new WATER();
1378 >      break;
1379  
1380 <  case FF_EAM:
1381 <    the_ff = new EAM_FF();
1382 <    break;
1383 <
1384 <  default:
1114 <    sprintf( painCave.errMsg,
1115 <       "SimSetup Error. Unrecognized force field in case statement.\n");
1116 <    painCave.isFatal = 1;
1117 <    simError();
1380 >    default:
1381 >      sprintf(painCave.errMsg,
1382 >              "SimSetup Error. Unrecognized force field in case statement.\n");
1383 >      painCave.isFatal = 1;
1384 >      simError();
1385    }
1386  
1387   #ifdef IS_MPI
1388 <  strcpy( checkPointMsg, "ForceField creation successful" );
1388 >  strcpy(checkPointMsg, "ForceField creation successful");
1389    MPIcheckPoint();
1390   #endif // is_mpi
1124 –
1391   }
1392  
1393  
1394 < void SimSetup::compList( void ){
1129 <
1394 > void SimSetup::compList(void){
1395    int i;
1396    char* id;
1397    LinkedMolStamp* headStamp = new LinkedMolStamp();
1398    LinkedMolStamp* currentStamp = NULL;
1399 <  comp_stamps = new MoleculeStamp*[n_components];
1399 >  comp_stamps = new MoleculeStamp * [n_components];
1400 >  bool haveCutoffGroups;
1401 >
1402 >  haveCutoffGroups = false;
1403    
1404    // make an array of molecule stamps that match the components used.
1405    // also extract the used stamps out into a separate linked list
1406 <  
1407 <  for(i=0; i<nInfo; i++){
1406 >
1407 >  for (i = 0; i < nInfo; i++){
1408      info[i].nComponents = n_components;
1409      info[i].componentsNmol = components_nmol;
1410      info[i].compStamps = comp_stamps;
1411      info[i].headStamp = headStamp;
1412    }
1145 –  
1413  
1147 –  for( i=0; i<n_components; i++ ){
1414  
1415 +  for (i = 0; i < n_components; i++){
1416      id = the_components[i]->getType();
1417      comp_stamps[i] = NULL;
1418 <    
1418 >
1419      // check to make sure the component isn't already in the list
1420  
1421 <    comp_stamps[i] = headStamp->match( id );
1422 <    if( comp_stamps[i] == NULL ){
1156 <      
1421 >    comp_stamps[i] = headStamp->match(id);
1422 >    if (comp_stamps[i] == NULL){
1423        // extract the component from the list;
1424 <      
1425 <      currentStamp = stamps->extractMolStamp( id );
1426 <      if( currentStamp == NULL ){
1427 <  sprintf( painCave.errMsg,
1428 <     "SimSetup error: Component \"%s\" was not found in the "
1429 <     "list of declared molecules\n",
1430 <     id );
1431 <  painCave.isFatal = 1;
1432 <  simError();
1424 >
1425 >      currentStamp = stamps->extractMolStamp(id);
1426 >      if (currentStamp == NULL){
1427 >        sprintf(painCave.errMsg,
1428 >                "SimSetup error: Component \"%s\" was not found in the "
1429 >                "list of declared molecules\n",
1430 >                id);
1431 >        painCave.isFatal = 1;
1432 >        simError();
1433        }
1434 <      
1435 <      headStamp->add( currentStamp );
1436 <      comp_stamps[i] = headStamp->match( id );
1434 >
1435 >      headStamp->add(currentStamp);
1436 >      comp_stamps[i] = headStamp->match(id);
1437      }
1438 +
1439 +    if(comp_stamps[i]->getNCutoffGroups() > 0)
1440 +      haveCutoffGroups = true;    
1441    }
1442 +    
1443 +  for (i = 0; i < nInfo; i++)
1444 +    info[i].haveCutoffGroups = haveCutoffGroups;
1445  
1446   #ifdef IS_MPI
1447 <  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
1447 >  strcpy(checkPointMsg, "Component stamps successfully extracted\n");
1448    MPIcheckPoint();
1449   #endif // is_mpi
1450 + }
1451  
1452 + void SimSetup::calcSysValues(void){
1453 +  int i;
1454  
1455 < }
1455 >  int* molMembershipArray;
1456  
1182 – void SimSetup::calcSysValues( void ){
1183 –  int i, j, k;
1184 –  
1185 –  int *molMembershipArray;
1186 –  
1457    tot_atoms = 0;
1458    tot_bonds = 0;
1459    tot_bends = 0;
1460    tot_torsions = 0;
1461 <  for( i=0; i<n_components; i++ ){
1462 <    
1463 <    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
1464 <    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
1465 <    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
1461 >  tot_rigid = 0;
1462 >  for (i = 0; i < n_components; i++){
1463 >    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
1464 >    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
1465 >    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
1466      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
1467 +    tot_rigid += components_nmol[i] * comp_stamps[i]->getNRigidBodies();
1468    }
1469    
1470    tot_SRI = tot_bonds + tot_bends + tot_torsions;
1471    molMembershipArray = new int[tot_atoms];
1472 <  
1473 <  for(i=0; i<nInfo; i++){
1472 >
1473 >  for (i = 0; i < nInfo; i++){
1474      info[i].n_atoms = tot_atoms;
1475      info[i].n_bonds = tot_bonds;
1476      info[i].n_bends = tot_bends;
1477      info[i].n_torsions = tot_torsions;
1478      info[i].n_SRI = tot_SRI;
1479      info[i].n_mol = tot_nmol;
1480 <    
1480 >
1481      info[i].molMembershipArray = molMembershipArray;
1482 <  }
1482 >  }
1483   }
1484  
1485   #ifdef IS_MPI
1486  
1487 < void SimSetup::mpiMolDivide( void ){
1217 <  
1487 > void SimSetup::mpiMolDivide(void){
1488    int i, j, k;
1489    int localMol, allMol;
1490    int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
1491 +  int local_rigid;
1492 +  vector<int> globalMolIndex;
1493  
1494 <  mpiSim = new mpiSimulation( info );
1223 <  
1224 <  globalIndex = mpiSim->divideLabor();
1494 >  mpiSim = new mpiSimulation(info);
1495  
1496 +  mpiSim->divideLabor();
1497 +  globalAtomIndex = mpiSim->getGlobalAtomIndex();
1498 +  //globalMolIndex = mpiSim->getGlobalMolIndex();
1499 +
1500    // set up the local variables
1501 <  
1501 >
1502    mol2proc = mpiSim->getMolToProcMap();
1503    molCompType = mpiSim->getMolComponentType();
1504 <  
1504 >
1505    allMol = 0;
1506    localMol = 0;
1507    local_atoms = 0;
1508    local_bonds = 0;
1509    local_bends = 0;
1510    local_torsions = 0;
1511 <  globalAtomIndex = 0;
1511 >  local_rigid = 0;
1512 >  globalAtomCounter = 0;
1513  
1514 <
1515 <  for( i=0; i<n_components; i++ ){
1516 <
1517 <    for( j=0; j<components_nmol[i]; j++ ){
1518 <      
1519 <      if( mol2proc[allMol] == worldRank ){
1520 <  
1521 <  local_atoms +=    comp_stamps[i]->getNAtoms();
1522 <  local_bonds +=    comp_stamps[i]->getNBonds();
1248 <  local_bends +=    comp_stamps[i]->getNBends();
1249 <  local_torsions += comp_stamps[i]->getNTorsions();
1250 <  localMol++;
1514 >  for (i = 0; i < n_components; i++){
1515 >    for (j = 0; j < components_nmol[i]; j++){
1516 >      if (mol2proc[allMol] == worldRank){
1517 >        local_atoms += comp_stamps[i]->getNAtoms();
1518 >        local_bonds += comp_stamps[i]->getNBonds();
1519 >        local_bends += comp_stamps[i]->getNBends();
1520 >        local_torsions += comp_stamps[i]->getNTorsions();
1521 >        local_rigid += comp_stamps[i]->getNRigidBodies();
1522 >        localMol++;
1523        }      
1524 <      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1525 <        info[0].molMembershipArray[globalAtomIndex] = allMol;
1526 <        globalAtomIndex++;
1524 >      for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1525 >        info[0].molMembershipArray[globalAtomCounter] = allMol;
1526 >        globalAtomCounter++;
1527        }
1528  
1529 <      allMol++;      
1529 >      allMol++;
1530      }
1531    }
1532    local_SRI = local_bonds + local_bends + local_torsions;
1533 <  
1533 >
1534    info[0].n_atoms = mpiSim->getMyNlocal();  
1535    
1536 <  if( local_atoms != info[0].n_atoms ){
1537 <    sprintf( painCave.errMsg,
1538 <       "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
1539 <       " localAtom (%d) are not equal.\n",
1540 <       info[0].n_atoms,
1541 <       local_atoms );
1536 >
1537 >  if (local_atoms != info[0].n_atoms){
1538 >    sprintf(painCave.errMsg,
1539 >            "SimSetup error: mpiSim's localAtom (%d) and SimSetup's\n"
1540 >            "\tlocalAtom (%d) are not equal.\n",
1541 >            info[0].n_atoms, local_atoms);
1542      painCave.isFatal = 1;
1543      simError();
1544    }
# Line 1277 | Line 1549 | void SimSetup::mpiMolDivide( void ){
1549    info[0].n_SRI = local_SRI;
1550    info[0].n_mol = localMol;
1551  
1552 <  strcpy( checkPointMsg, "Passed nlocal consistency check." );
1552 >  strcpy(checkPointMsg, "Passed nlocal consistency check.");
1553    MPIcheckPoint();
1554   }
1555 <
1555 >
1556   #endif // is_mpi
1557  
1558  
1559 < void SimSetup::makeSysArrays( void ){
1560 <  int i, j, k, l;
1559 > void SimSetup::makeSysArrays(void){
1560 >
1561 > #ifndef IS_MPI
1562 >  int k, j;
1563 > #endif // is_mpi
1564 >  int i, l;
1565  
1566    Atom** the_atoms;
1567    Molecule* the_molecules;
1292 –  Exclude** the_excludes;
1568  
1569 <  
1295 <  for(l=0; l<nInfo; l++){
1296 <    
1569 >  for (l = 0; l < nInfo; l++){
1570      // create the atom and short range interaction arrays
1571 <    
1572 <    the_atoms = new Atom*[info[l].n_atoms];
1571 >
1572 >    the_atoms = new Atom * [info[l].n_atoms];
1573      the_molecules = new Molecule[info[l].n_mol];
1574      int molIndex;
1575  
1576      // initialize the molecule's stampID's
1577 <    
1577 >
1578   #ifdef IS_MPI
1579 <    
1580 <    
1579 >
1580 >
1581      molIndex = 0;
1582 <    for(i=0; i<mpiSim->getTotNmol(); i++){
1583 <    
1584 <      if(mol2proc[i] == worldRank ){
1585 <  the_molecules[molIndex].setStampID( molCompType[i] );
1586 <  the_molecules[molIndex].setMyIndex( molIndex );
1587 <  the_molecules[molIndex].setGlobalIndex( i );
1315 <  molIndex++;
1582 >    for (i = 0; i < mpiSim->getTotNmol(); i++){
1583 >      if (mol2proc[i] == worldRank){
1584 >        the_molecules[molIndex].setStampID(molCompType[i]);
1585 >        the_molecules[molIndex].setMyIndex(molIndex);
1586 >        the_molecules[molIndex].setGlobalIndex(i);
1587 >        molIndex++;
1588        }
1589      }
1590 <    
1590 >
1591   #else // is_mpi
1592 <    
1592 >
1593      molIndex = 0;
1594 <    globalAtomIndex = 0;
1595 <    for(i=0; i<n_components; i++){
1596 <      for(j=0; j<components_nmol[i]; j++ ){
1597 <  the_molecules[molIndex].setStampID( i );
1598 <  the_molecules[molIndex].setMyIndex( molIndex );
1599 <  the_molecules[molIndex].setGlobalIndex( molIndex );
1600 <  for (k = 0; k < comp_stamps[i]->getNAtoms(); k++) {
1601 <    info[l].molMembershipArray[globalAtomIndex] = molIndex;
1602 <    globalAtomIndex++;
1603 <  }
1604 <  molIndex++;
1594 >    globalAtomCounter = 0;
1595 >    for (i = 0; i < n_components; i++){
1596 >      for (j = 0; j < components_nmol[i]; j++){
1597 >        the_molecules[molIndex].setStampID(i);
1598 >        the_molecules[molIndex].setMyIndex(molIndex);
1599 >        the_molecules[molIndex].setGlobalIndex(molIndex);
1600 >        for (k = 0; k < comp_stamps[i]->getNAtoms(); k++){
1601 >          info[l].molMembershipArray[globalAtomCounter] = molIndex;
1602 >          globalAtomCounter++;
1603 >        }
1604 >        molIndex++;
1605        }
1606      }
1335 –    
1336 –    
1337 – #endif // is_mpi
1607  
1608  
1609 <    if( info[l].n_SRI ){
1341 <    
1342 <      Exclude::createArray(info[l].n_SRI);
1343 <      the_excludes = new Exclude*[info[l].n_SRI];
1344 <      for( int ex=0; ex<info[l].n_SRI; ex++){
1345 <  the_excludes[ex] = new Exclude(ex);
1346 <      }
1347 <      info[l].globalExcludes = new int;
1348 <      info[l].n_exclude = info[l].n_SRI;
1349 <    }
1350 <    else{
1351 <    
1352 <      Exclude::createArray( 1 );
1353 <      the_excludes = new Exclude*;
1354 <      the_excludes[0] = new Exclude(0);
1355 <      the_excludes[0]->setPair( 0,0 );
1356 <      info[l].globalExcludes = new int;
1357 <      info[l].globalExcludes[0] = 0;
1358 <      info[l].n_exclude = 0;
1359 <    }
1609 > #endif // is_mpi
1610  
1611 +    info[l].globalExcludes = new int;
1612 +    info[l].globalExcludes[0] = 0;
1613 +    
1614      // set the arrays into the SimInfo object
1615  
1616      info[l].atoms = the_atoms;
1617      info[l].molecules = the_molecules;
1618      info[l].nGlobalExcludes = 0;
1366 –    info[l].excludes = the_excludes;
1367 –
1368 –    the_ff->setSimInfo( info );
1619      
1620 +    the_ff->setSimInfo(info);
1621    }
1622   }
1623  
1624 < void SimSetup::makeIntegrator( void ){
1374 <
1624 > void SimSetup::makeIntegrator(void){
1625    int k;
1626  
1627 <  NVT<RealIntegrator>*  myNVT = NULL;
1628 <  NPTi<RealIntegrator>* myNPTi = NULL;
1629 <  NPTf<RealIntegrator>* myNPTf = NULL;
1630 <  NPTim<RealIntegrator>* myNPTim = NULL;
1631 <  NPTfm<RealIntegrator>* myNPTfm = NULL;
1632 <        
1633 <  for(k=0; k<nInfo; k++){
1634 <    
1635 <    switch( ensembleCase ){
1636 <      
1637 <    case NVE_ENS:
1638 <      if (globals->haveZconstraints()){
1639 <        setupZConstraint(info[k]);
1640 <        new ZConstraint<NVE<RealIntegrator> >( &(info[k]), the_ff );
1641 <     }
1627 >  NVE<RealIntegrator>* myNVE = NULL;
1628 >  NVT<RealIntegrator>* myNVT = NULL;
1629 >  NPTi<NPT<RealIntegrator> >* myNPTi = NULL;
1630 >  NPTf<NPT<RealIntegrator> >* myNPTf = NULL;
1631 >  NPTxyz<NPT<RealIntegrator> >* myNPTxyz = NULL;
1632 >  
1633 >  for (k = 0; k < nInfo; k++){
1634 >    switch (ensembleCase){
1635 >      case NVE_ENS:
1636 >        if (globals->haveZconstraints()){
1637 >          setupZConstraint(info[k]);
1638 >          myNVE = new ZConstraint<NVE<RealIntegrator> >(&(info[k]), the_ff);
1639 >        }
1640 >        else{
1641 >          myNVE = new NVE<RealIntegrator>(&(info[k]), the_ff);
1642 >        }
1643 >        
1644 >        info->the_integrator = myNVE;
1645 >        break;
1646  
1647 <     else
1648 <        new NVE<RealIntegrator>( &(info[k]), the_ff );
1649 <      break;
1650 <      
1651 <    case NVT_ENS:
1652 <      if (globals->haveZconstraints()){
1653 <        setupZConstraint(info[k]);
1400 <        myNVT = new ZConstraint<NVT<RealIntegrator> >( &(info[k]), the_ff );
1401 <      }
1402 <      else
1403 <        myNVT = new NVT<RealIntegrator>( &(info[k]), the_ff );
1647 >      case NVT_ENS:
1648 >        if (globals->haveZconstraints()){
1649 >          setupZConstraint(info[k]);
1650 >          myNVT = new ZConstraint<NVT<RealIntegrator> >(&(info[k]), the_ff);
1651 >        }
1652 >        else
1653 >          myNVT = new NVT<RealIntegrator>(&(info[k]), the_ff);
1654  
1655          myNVT->setTargetTemp(globals->getTargetTemp());
1656 <      
1657 <        if (globals->haveTauThermostat())
1656 >
1657 >        if (globals->haveTauThermostat())
1658            myNVT->setTauThermostat(globals->getTauThermostat());
1659 <      
1660 <        else {
1661 <          sprintf( painCave.errMsg,
1662 <                    "SimSetup error: If you use the NVT\n"
1413 <                    "    ensemble, you must set tauThermostat.\n");
1659 >        else{
1660 >          sprintf(painCave.errMsg,
1661 >                  "SimSetup error: If you use the NVT\n"
1662 >                  "\tensemble, you must set tauThermostat.\n");
1663            painCave.isFatal = 1;
1664            simError();
1665          }
1666 +
1667 +        info->the_integrator = myNVT;
1668          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 );
1669  
1670 <      myNPTi->setTargetTemp( globals->getTargetTemp() );
1671 <          
1672 <      if (globals->haveTargetPressure())
1673 <        myNPTi->setTargetPressure(globals->getTargetPressure());
1674 <      else {
1675 <         sprintf( painCave.errMsg,
1676 <                   "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 );
1670 >      case NPTi_ENS:
1671 >        if (globals->haveZconstraints()){
1672 >          setupZConstraint(info[k]);
1673 >          myNPTi = new ZConstraint<NPTi<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1674 >        }
1675 >        else
1676 >          myNPTi = new NPTi<NPT<RealIntegrator> >(&(info[k]), the_ff);
1677  
1678 <      myNPTf->setTargetTemp( globals->getTargetTemp());
1679 <          
1680 <      if (globals->haveTargetPressure())
1681 <        myNPTf->setTargetPressure(globals->getTargetPressure());
1682 <      else {
1683 <        sprintf( painCave.errMsg,
1678 >        myNPTi->setTargetTemp(globals->getTargetTemp());
1679 >
1680 >        if (globals->haveTargetPressure())
1681 >          myNPTi->setTargetPressure(globals->getTargetPressure());
1682 >        else{
1683 >          sprintf(painCave.errMsg,
1684                    "SimSetup error: If you use a constant pressure\n"
1685 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1686 <        painCave.isFatal = 1;
1687 <        simError();
1688 <      }    
1689 <          
1690 <      if( globals->haveTauThermostat() )
1691 <        myNPTf->setTauThermostat( globals->getTauThermostat() );
1692 <      else{
1693 <        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,
1685 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1686 >          painCave.isFatal = 1;
1687 >          simError();
1688 >        }
1689 >
1690 >        if (globals->haveTauThermostat())
1691 >          myNPTi->setTauThermostat(globals->getTauThermostat());
1692 >        else{
1693 >          sprintf(painCave.errMsg,
1694                    "SimSetup error: If you use an NPT\n"
1695 <                  "    ensemble, you must set tauBarostat.\n");
1696 <        painCave.isFatal = 1;
1697 <        simError();
1698 <      }
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 );
1695 >                  "\tensemble, you must set tauThermostat.\n");
1696 >          painCave.isFatal = 1;
1697 >          simError();
1698 >        }
1699  
1700 <        myNPTim->setTargetTemp( globals->getTargetTemp());
1701 <          
1702 <      if (globals->haveTargetPressure())
1703 <        myNPTim->setTargetPressure(globals->getTargetPressure());
1513 <      else {
1514 <        sprintf( painCave.errMsg,
1515 <                  "SimSetup error: If you use a constant pressure\n"
1516 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1517 <        painCave.isFatal = 1;
1518 <        simError();
1519 <      }
1520 <          
1521 <      if( globals->haveTauThermostat() )
1522 <        myNPTim->setTauThermostat( globals->getTauThermostat() );
1523 <      else{
1524 <        sprintf( painCave.errMsg,
1700 >        if (globals->haveTauBarostat())
1701 >          myNPTi->setTauBarostat(globals->getTauBarostat());
1702 >        else{
1703 >          sprintf(painCave.errMsg,
1704                    "SimSetup error: If you use an NPT\n"
1705 <                  "    ensemble, you must set tauThermostat.\n");
1706 <        painCave.isFatal = 1;
1707 <        simError();
1708 <      }
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 );
1705 >                  "\tensemble, you must set tauBarostat.\n");
1706 >          painCave.isFatal = 1;
1707 >          simError();
1708 >        }
1709  
1710 <      myNPTfm->setTargetTemp( globals->getTargetTemp());
1710 >        info->the_integrator = myNPTi;
1711 >        break;
1712  
1713 <      if (globals->haveTargetPressure())
1714 <        myNPTfm->setTargetPressure(globals->getTargetPressure());
1715 <      else {
1716 <        sprintf( painCave.errMsg,
1713 >      case NPTf_ENS:
1714 >        if (globals->haveZconstraints()){
1715 >          setupZConstraint(info[k]);
1716 >          myNPTf = new ZConstraint<NPTf<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1717 >        }
1718 >        else
1719 >          myNPTf = new NPTf<NPT <RealIntegrator> >(&(info[k]), the_ff);
1720 >
1721 >        myNPTf->setTargetTemp(globals->getTargetTemp());
1722 >
1723 >        if (globals->haveTargetPressure())
1724 >          myNPTf->setTargetPressure(globals->getTargetPressure());
1725 >        else{
1726 >          sprintf(painCave.errMsg,
1727                    "SimSetup error: If you use a constant pressure\n"
1728 <                  "    ensemble, you must set targetPressure in the BASS file.\n");
1729 <        painCave.isFatal = 1;
1730 <        simError();
1731 <      }
1728 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1729 >          painCave.isFatal = 1;
1730 >          simError();
1731 >        }    
1732  
1733 <      if( globals->haveTauThermostat() )
1734 <        myNPTfm->setTauThermostat( globals->getTauThermostat() );
1735 <      else{
1736 <        sprintf( painCave.errMsg,
1733 >        if (globals->haveTauThermostat())
1734 >          myNPTf->setTauThermostat(globals->getTauThermostat());
1735 >
1736 >        else{
1737 >          sprintf(painCave.errMsg,
1738                    "SimSetup error: If you use an NPT\n"
1739 <                  "    ensemble, you must set tauThermostat.\n");
1740 <        painCave.isFatal = 1;
1741 <        simError();
1742 <      }
1739 >                  "\tensemble, you must set tauThermostat.\n");
1740 >          painCave.isFatal = 1;
1741 >          simError();
1742 >        }
1743  
1744 <      if( globals->haveTauBarostat() )
1745 <        myNPTfm->setTauBarostat( globals->getTauBarostat() );
1746 <      else{
1747 <        sprintf( painCave.errMsg,
1744 >        if (globals->haveTauBarostat())
1745 >          myNPTf->setTauBarostat(globals->getTauBarostat());
1746 >
1747 >        else{
1748 >          sprintf(painCave.errMsg,
1749                    "SimSetup error: If you use an NPT\n"
1750 <                  "    ensemble, you must set tauBarostat.\n");
1750 >                  "\tensemble, you must set tauBarostat.\n");
1751 >          painCave.isFatal = 1;
1752 >          simError();
1753 >        }
1754 >
1755 >        info->the_integrator = myNPTf;
1756 >        break;
1757 >
1758 >      case NPTxyz_ENS:
1759 >        if (globals->haveZconstraints()){
1760 >          setupZConstraint(info[k]);
1761 >          myNPTxyz = new ZConstraint<NPTxyz<NPT <RealIntegrator> > >(&(info[k]), the_ff);
1762 >        }
1763 >        else
1764 >          myNPTxyz = new NPTxyz<NPT <RealIntegrator> >(&(info[k]), the_ff);
1765 >
1766 >        myNPTxyz->setTargetTemp(globals->getTargetTemp());
1767 >
1768 >        if (globals->haveTargetPressure())
1769 >          myNPTxyz->setTargetPressure(globals->getTargetPressure());
1770 >        else{
1771 >          sprintf(painCave.errMsg,
1772 >                  "SimSetup error: If you use a constant pressure\n"
1773 >                  "\tensemble, you must set targetPressure in the BASS file.\n");
1774 >          painCave.isFatal = 1;
1775 >          simError();
1776 >        }    
1777 >
1778 >        if (globals->haveTauThermostat())
1779 >          myNPTxyz->setTauThermostat(globals->getTauThermostat());
1780 >        else{
1781 >          sprintf(painCave.errMsg,
1782 >                  "SimSetup error: If you use an NPT\n"
1783 >                  "\tensemble, you must set tauThermostat.\n");
1784 >          painCave.isFatal = 1;
1785 >          simError();
1786 >        }
1787 >
1788 >        if (globals->haveTauBarostat())
1789 >          myNPTxyz->setTauBarostat(globals->getTauBarostat());
1790 >        else{
1791 >          sprintf(painCave.errMsg,
1792 >                  "SimSetup error: If you use an NPT\n"
1793 >                  "\tensemble, you must set tauBarostat.\n");
1794 >          painCave.isFatal = 1;
1795 >          simError();
1796 >        }
1797 >
1798 >        info->the_integrator = myNPTxyz;
1799 >        break;
1800 >
1801 >      default:
1802 >        sprintf(painCave.errMsg,
1803 >                "SimSetup Error. Unrecognized ensemble in case statement.\n");
1804          painCave.isFatal = 1;
1805          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();
1806      }
1807    }
1808   }
1809  
1810 < void SimSetup::initFortran( void ){
1593 <
1810 > void SimSetup::initFortran(void){
1811    info[0].refreshSim();
1812 <  
1813 <  if( !strcmp( info[0].mixingRule, "standard") ){
1814 <    the_ff->initForceField( LB_MIXING_RULE );
1812 >
1813 >  if (!strcmp(info[0].mixingRule, "standard")){
1814 >    the_ff->initForceField(LB_MIXING_RULE);
1815    }
1816 <  else if( !strcmp( info[0].mixingRule, "explicit") ){
1817 <    the_ff->initForceField( EXPLICIT_MIXING_RULE );
1816 >  else if (!strcmp(info[0].mixingRule, "explicit")){
1817 >    the_ff->initForceField(EXPLICIT_MIXING_RULE);
1818    }
1819    else{
1820 <    sprintf( painCave.errMsg,
1821 <       "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1605 <       info[0].mixingRule );
1820 >    sprintf(painCave.errMsg, "SimSetup Error: unknown mixing rule -> \"%s\"\n",
1821 >            info[0].mixingRule);
1822      painCave.isFatal = 1;
1823      simError();
1824    }
1825  
1826  
1827   #ifdef IS_MPI
1828 <  strcpy( checkPointMsg,
1613 <    "Successfully intialized the mixingRule for Fortran." );
1828 >  strcpy(checkPointMsg, "Successfully intialized the mixingRule for Fortran.");
1829    MPIcheckPoint();
1830   #endif // is_mpi
1616 –
1831   }
1832  
1833 < void SimSetup::setupZConstraint(SimInfo& theInfo)
1620 < {
1833 > void SimSetup::setupZConstraint(SimInfo& theInfo){
1834    int nZConstraints;
1835    ZconStamp** zconStamp;
1836  
1837 <  if(globals->haveZconstraintTime()){  
1625 <    
1837 >  if (globals->haveZconstraintTime()){
1838      //add sample time of z-constraint  into SimInfo's property list                    
1839      DoubleData* zconsTimeProp = new DoubleData();
1840      zconsTimeProp->setID(ZCONSTIME_ID);
# Line 1630 | Line 1842 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1842      theInfo.addProperty(zconsTimeProp);
1843    }
1844    else{
1845 <    sprintf( painCave.errMsg,
1846 <       "ZConstraint error: If you use an ZConstraint\n"
1847 <       " , you must set sample time.\n");
1845 >    sprintf(painCave.errMsg,
1846 >            "ZConstraint error: If you use a ZConstraint,\n"
1847 >            "\tyou must set zconsTime.\n");
1848      painCave.isFatal = 1;
1849 <    simError();      
1849 >    simError();
1850    }
1851  
1852    //push zconsTol into siminfo, if user does not specify
1853    //value for zconsTol, a default value will be used
1854    DoubleData* zconsTol = new DoubleData();
1855    zconsTol->setID(ZCONSTOL_ID);
1856 <  if(globals->haveZconsTol()){
1856 >  if (globals->haveZconsTol()){
1857      zconsTol->setData(globals->getZconsTol());
1858    }
1859    else{
1860 <  double defaultZConsTol = 0.01;
1861 <    sprintf( painCave.errMsg,
1862 <       "ZConstraint Waring: Tolerance for z-constraint methodl is not specified\n"
1863 <       " , default value %f is used.\n", defaultZConsTol);
1860 >    double defaultZConsTol = 0.01;
1861 >    sprintf(painCave.errMsg,
1862 >            "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n"
1863 >            "\tOOPSE will use a default value of %f.\n"
1864 >            "\tTo set the tolerance, use the zconsTol variable.\n",
1865 >            defaultZConsTol);
1866      painCave.isFatal = 0;
1867      simError();      
1868  
# Line 1656 | Line 1870 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1870    }
1871    theInfo.addProperty(zconsTol);
1872  
1873 <  //set Force Substraction Policy
1874 <  StringData* zconsForcePolicy =  new StringData();
1873 >  //set Force Subtraction Policy
1874 >  StringData* zconsForcePolicy = new StringData();
1875    zconsForcePolicy->setID(ZCONSFORCEPOLICY_ID);
1876 <  
1877 <  if(globals->haveZconsForcePolicy()){
1876 >
1877 >  if (globals->haveZconsForcePolicy()){
1878      zconsForcePolicy->setData(globals->getZconsForcePolicy());
1879 <  }  
1879 >  }
1880    else{
1881 <     sprintf( painCave.errMsg,
1882 <             "ZConstraint Warning: User does not set force substraction policy, "
1883 <             "average force substraction policy is used\n");
1884 <     painCave.isFatal = 0;
1885 <     simError();
1886 <     zconsForcePolicy->setData("BYNUMBER");
1881 >    sprintf(painCave.errMsg,
1882 >            "ZConstraint Warning: No force subtraction policy was set.\n"
1883 >            "\tOOPSE will use PolicyByMass.\n"
1884 >            "\tTo set the policy, use the zconsForcePolicy variable.\n");
1885 >    painCave.isFatal = 0;
1886 >    simError();
1887 >    zconsForcePolicy->setData("BYMASS");
1888    }
1889 <
1890 < theInfo.addProperty(zconsForcePolicy);
1891 <
1889 >
1890 >  theInfo.addProperty(zconsForcePolicy);
1891 >
1892 >  //set zcons gap
1893 >  DoubleData* zconsGap = new DoubleData();
1894 >  zconsGap->setID(ZCONSGAP_ID);
1895 >
1896 >  if (globals->haveZConsGap()){
1897 >    zconsGap->setData(globals->getZconsGap());
1898 >    theInfo.addProperty(zconsGap);  
1899 >  }
1900 >
1901 >  //set zcons fixtime
1902 >  DoubleData* zconsFixtime = new DoubleData();
1903 >  zconsFixtime->setID(ZCONSFIXTIME_ID);
1904 >
1905 >  if (globals->haveZConsFixTime()){
1906 >    zconsFixtime->setData(globals->getZconsFixtime());
1907 >    theInfo.addProperty(zconsFixtime);  
1908 >  }
1909 >
1910 >  //set zconsUsingSMD
1911 >  IntData* zconsUsingSMD = new IntData();
1912 >  zconsUsingSMD->setID(ZCONSUSINGSMD_ID);
1913 >
1914 >  if (globals->haveZConsUsingSMD()){
1915 >    zconsUsingSMD->setData(globals->getZconsUsingSMD());
1916 >    theInfo.addProperty(zconsUsingSMD);  
1917 >  }
1918 >
1919    //Determine the name of ouput file and add it into SimInfo's property list
1920    //Be careful, do not use inFileName, since it is a pointer which
1921    //point to a string at master node, and slave nodes do not contain that string
1922 <  
1922 >
1923    string zconsOutput(theInfo.finalName);
1924 <  
1924 >
1925    zconsOutput = zconsOutput.substr(0, zconsOutput.rfind(".")) + ".fz";
1926 <  
1926 >
1927    StringData* zconsFilename = new StringData();
1928    zconsFilename->setID(ZCONSFILENAME_ID);
1929    zconsFilename->setData(zconsOutput);
1930 <  
1930 >
1931    theInfo.addProperty(zconsFilename);
1932 <  
1932 >
1933    //setup index, pos and other parameters of z-constraint molecules
1934    nZConstraints = globals->getNzConstraints();
1935    theInfo.nZconstraints = nZConstraints;
# Line 1697 | Line 1939 | void SimSetup::setupZConstraint(SimInfo& theInfo)
1939  
1940    ZConsParaData* zconsParaData = new ZConsParaData();
1941    zconsParaData->setID(ZCONSPARADATA_ID);
1942 <
1943 <  for(int i = 0; i < nZConstraints; i++){
1942 >
1943 >  for (int i = 0; i < nZConstraints; i++){
1944      tempParaItem.havingZPos = zconStamp[i]->haveZpos();
1945      tempParaItem.zPos = zconStamp[i]->getZpos();
1946      tempParaItem.zconsIndex = zconStamp[i]->getMolIndex();
1947      tempParaItem.kRatio = zconStamp[i]->getKratio();
1948 <
1948 >    tempParaItem.havingCantVel = zconStamp[i]->haveCantVel();
1949 >    tempParaItem.cantVel = zconStamp[i]->getCantVel();    
1950      zconsParaData->addItem(tempParaItem);
1951    }
1952  
1953 +  //check the uniqueness of index  
1954 +  if(!zconsParaData->isIndexUnique()){
1955 +    sprintf(painCave.errMsg,
1956 +            "ZConstraint Error: molIndex is not unique!\n");
1957 +    painCave.isFatal = 1;
1958 +    simError();
1959 +  }
1960 +
1961    //sort the parameters by index of molecules
1962    zconsParaData->sortByIndex();
1963    
1964    //push data into siminfo, therefore, we can retrieve later
1965    theInfo.addProperty(zconsParaData);
1715 –      
1966   }
1967 +
1968 + void SimSetup::makeMinimizer(){
1969 +
1970 +  OOPSEMinimizer* myOOPSEMinimizer;
1971 +  MinimizerParameterSet* param;
1972 +  char minimizerName[100];
1973 +  
1974 +  for (int i = 0; i < nInfo; i++){
1975 +    
1976 +    //prepare parameter set for minimizer
1977 +    param = new MinimizerParameterSet();
1978 +    param->setDefaultParameter();
1979 +
1980 +    if (globals->haveMinimizer()){
1981 +      param->setFTol(globals->getMinFTol());
1982 +    }
1983 +
1984 +    if (globals->haveMinGTol()){
1985 +      param->setGTol(globals->getMinGTol());
1986 +    }
1987 +
1988 +    if (globals->haveMinMaxIter()){
1989 +      param->setMaxIteration(globals->getMinMaxIter());
1990 +    }
1991 +
1992 +    if (globals->haveMinWriteFrq()){
1993 +      param->setMaxIteration(globals->getMinMaxIter());
1994 +    }
1995 +
1996 +    if (globals->haveMinWriteFrq()){
1997 +      param->setWriteFrq(globals->getMinWriteFrq());
1998 +    }
1999 +    
2000 +    if (globals->haveMinStepSize()){
2001 +      param->setStepSize(globals->getMinStepSize());
2002 +    }
2003 +
2004 +    if (globals->haveMinLSMaxIter()){
2005 +      param->setLineSearchMaxIteration(globals->getMinLSMaxIter());
2006 +    }    
2007 +
2008 +    if (globals->haveMinLSTol()){
2009 +      param->setLineSearchTol(globals->getMinLSTol());
2010 +    }    
2011 +
2012 +    strcpy(minimizerName, globals->getMinimizer());
2013 +
2014 +    if (!strcasecmp(minimizerName, "CG")){
2015 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);
2016 +    }
2017 +    else if (!strcasecmp(minimizerName, "SD")){
2018 +    //myOOPSEMinimizer = MinimizerFactory.creatMinimizer("", &(info[i]), the_ff, param);
2019 +      myOOPSEMinimizer = new SDMinimizer(&(info[i]), the_ff, param);
2020 +    }
2021 +    else{
2022 +          sprintf(painCave.errMsg,
2023 +                  "SimSetup error: Unrecognized Minimizer, use Conjugate Gradient \n");
2024 +          painCave.isFatal = 0;
2025 +          simError();
2026 +
2027 +      myOOPSEMinimizer = new PRCGMinimizer(&(info[i]), the_ff, param);          
2028 +    }
2029 +     info[i].the_integrator = myOOPSEMinimizer;
2030 +
2031 +     //store the minimizer into simInfo
2032 +     info[i].the_minimizer = myOOPSEMinimizer;
2033 +     info[i].has_minimizer = true;
2034 +  }
2035 +
2036 + }

Diff Legend

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