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 614 by mmeineke, Tue Jul 15 17:57:04 2003 UTC vs.
Revision 708 by tim, Wed Aug 20 22:23:34 2003 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines