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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines