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 812 by mmeineke, Wed Oct 22 21:17:32 2003 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines