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 412 by mmeineke, Wed Mar 26 21:50:33 2003 UTC vs.
Revision 965 by gezelter, Mon Jan 19 21:17:39 2004 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines