ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/SimSetup.cpp
(Generate patch)

Comparing trunk/mdtools/interface_implementation/SimSetup.cpp (file contents):
Revision 113 by mmeineke, Mon Sep 23 15:12:56 2002 UTC vs.
Revision 205 by mmeineke, Wed Dec 11 20:39:41 2002 UTC

# Line 6 | Line 6
6   #include "parse_me.h"
7   #include "LRI.hpp"
8   #include "Integrator.hpp"
9 + #include "simError.h"
10  
11 + #ifdef IS_MPI
12 + #include "mpiBASS.h"
13 + #include "bassDiag.hpp"
14 + #endif
15 +
16   SimSetup::SimSetup(){
17    stamps = new MakeStamps();
18    globals = new Globals();
19 +  
20 + #ifdef IS_MPI
21 +  strcpy( checkPointMsg, "SimSetup creation successful" );
22 +  MPIcheckPoint();
23 + #endif // IS_MPI
24   }
25  
26   SimSetup::~SimSetup(){
# Line 19 | Line 30 | void SimSetup::parseFile( char* fileName ){
30  
31   void SimSetup::parseFile( char* fileName ){
32  
33 <  inFileName = fileName;
34 <  set_interface_stamps( stamps, globals );
35 <  yacc_BASS( fileName );
33 > #ifdef IS_MPI
34 >  if( worldRank == 0 ){
35 > #endif // is_mpi
36 >    
37 >    inFileName = fileName;
38 >    set_interface_stamps( stamps, globals );
39 >    
40 > #ifdef IS_MPI
41 >    mpiEventInit();
42 > #endif
43 >
44 >    yacc_BASS( fileName );
45 >
46 > #ifdef IS_MPI
47 >    throwMPIEvent(NULL);
48 >  }
49 >  else receiveParse();
50 > #endif
51 >
52   }
53  
54 + #ifdef IS_MPI
55 + void SimSetup::receiveParse(void){
56 +
57 +    set_interface_stamps( stamps, globals );
58 +    mpiEventInit();
59 +    MPIcheckPoint();
60 +    mpiEventLoop();
61 +
62 + }
63 +
64 +
65 + void SimSetup::testMe(void){
66 +  bassDiag* dumpMe = new bassDiag(globals,stamps);
67 +  dumpMe->dumpStamps();
68 +  delete dumpMe;
69 + }
70 + #endif
71 +
72   void SimSetup::createSim( void ){
73  
74    MakeStamps *the_stamps;
75    Globals* the_globals;
76 <  int i;
76 >  int i, j;
77  
78    // get the stamps and globals;
79    the_stamps = stamps;
# Line 43 | Line 88 | void SimSetup::createSim( void ){
88    n_components = the_globals->getNComponents();
89    strcpy( force_field, the_globals->getForceField() );
90    strcpy( ensemble, the_globals->getEnsemble() );
91 <  
91 >
92    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
93    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
94    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
95    else{
96 <    std::cerr<< "SimSetup Error. Unrecognized force field -> "
97 <             << force_field << "\n";
98 <    exit(8);
96 >    sprintf( painCave.errMsg,
97 >             "SimSetup Error. Unrecognized force field -> %s\n",
98 >             force_field );
99 >    painCave.isFatal = 1;
100 >    simError();
101    }
102  
103 + #ifdef IS_MPI
104 +  strcpy( checkPointMsg, "ForceField creation successful" );
105 +  MPIcheckPoint();
106 + #endif // is_mpi
107 +
108    // get the components and calculate the tot_nMol and indvidual n_mol
109    the_components = the_globals->getComponents();
110    components_nmol = new int[n_components];
111    comp_stamps = new MoleculeStamp*[n_components];
112 <  
112 >
113    if( !the_globals->haveNMol() ){
114 <    // we don't have the total number of molecules, so we assume it is
114 >    // we don't have the total number of molecules, so we assume it is
115      // given in each component
116  
117      tot_nmol = 0;
118      for( i=0; i<n_components; i++ ){
119 <      
119 >
120        if( !the_components[i]->haveNMol() ){
121          // we have a problem
122 <        std::cerr << "SimSetup Error. No global NMol or component NMol"
123 <                  << " given. Cannot calculate the number of atoms.\n";
124 <        exit( 8 );
122 >        sprintf( painCave.errMsg,
123 >                 "SimSetup Error. No global NMol or component NMol"
124 >                 " given. Cannot calculate the number of atoms.\n" );
125 >        painCave.isFatal = 1;
126 >        simError();
127        }
128  
129        tot_nmol += the_components[i]->getNMol();
# Line 77 | Line 131 | void SimSetup::createSim( void ){
131      }
132    }
133    else{
134 <    std::cerr << "NOT A SUPPORTED FEATURE\n";
134 >    sprintf( painCave.errMsg,
135 >             "SimSetup error.\n"
136 >             "\tSorry, the ability to specify total"
137 >             " nMols and then give molfractions in the components\n"
138 >             "\tis not currently supported."
139 >             " Please give nMol in the components.\n" );
140 >    painCave.isFatal = 1;
141 >    simError();
142      
82 – //     tot_nmol = the_globals->getNMol();
143      
144 < //     //we have the total number of molecules, now we check for molfractions
145 < //     for( i=0; i<n_components; i++ ){
146 <      
147 < //       if( !the_components[i]->haveMolFraction() ){
148 <        
149 < //      if( !the_components[i]->haveNMol() ){
150 < //        //we have a problem
151 < //        std::cerr << "SimSetup error. Neither molFraction nor "
152 < //                  << " nMol was given in component
153 <
144 >    //     tot_nmol = the_globals->getNMol();
145 >    
146 >    //   //we have the total number of molecules, now we check for molfractions
147 >    //     for( i=0; i<n_components; i++ ){
148 >    
149 >    //       if( !the_components[i]->haveMolFraction() ){
150 >    
151 >    //  if( !the_components[i]->haveNMol() ){
152 >    //    //we have a problem
153 >    //    std::cerr << "SimSetup error. Neither molFraction nor "
154 >    //              << " nMol was given in component
155 >    
156    }
157  
158 + #ifdef IS_MPI
159 +  strcpy( checkPointMsg, "Have the number of components" );
160 +  MPIcheckPoint();
161 + #endif // is_mpi
162 +
163    // make an array of molecule stamps that match the components used.
164 +  // also extract the used stamps out into a separate linked list
165  
166 +  simnfo->nComponents = n_components;
167 +  simnfo->componentsNmol = components_nmol;
168 +  simnfo->compStamps = comp_stamps;
169 +  simnfo->headStamp = new LinkedMolStamp();
170 +  
171 +  char* id;
172 +  LinkedMolStamp* headStamp = simnfo->headStamp;
173 +  LinkedMolStamp* currentStamp = NULL;
174    for( i=0; i<n_components; i++ ){
175  
176 <    comp_stamps[i] =
177 <      the_stamps->getMolecule( the_components[i]->getType() );
176 >    id = the_components[i]->getType();
177 >    comp_stamps[i] = NULL;
178 >    
179 >    // check to make sure the component isn't already in the list
180 >
181 >    comp_stamps[i] = headStamp->match( id );
182 >    if( comp_stamps[i] == NULL ){
183 >      
184 >      // extract the component from the list;
185 >      
186 >      currentStamp = the_stamps->extractMolStamp( id );
187 >      if( currentStamp == NULL ){
188 >        sprintf( painCave.errMsg,
189 >                 "SimSetup error: Component \"%s\" was not found in the "
190 >                 "list of declared molecules\n"
191 >                 id );
192 >        painCave.isFatal = 1;
193 >        simError();
194 >      }
195 >      
196 >      headStamp->add( currentStamp );
197 >      comp_stamps[i] = headStamp->match( id );
198 >    }
199    }
200  
201 + #ifdef IS_MPI
202 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
203 +  MPIcheckPoint();
204 + #endif // is_mpi
205    
206  
207 +
208 +
209    // caclulate the number of atoms, bonds, bends and torsions
210  
211    tot_atoms = 0;
# Line 111 | Line 214 | void SimSetup::createSim( void ){
214    tot_torsions = 0;
215    for( i=0; i<n_components; i++ ){
216      
217 <    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
218 <    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
219 <    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
217 >    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
218 >    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
219 >    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
220      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
221    }
222 <  
222 >
223    tot_SRI = tot_bonds + tot_bends + tot_torsions;
224 <  
224 >
225    simnfo->n_atoms = tot_atoms;
226    simnfo->n_bonds = tot_bonds;
227    simnfo->n_bends = tot_bends;
228    simnfo->n_torsions = tot_torsions;
229    simnfo->n_SRI = tot_SRI;
230 +  simnfo->n_mol = tot_nmol;
231  
128 –  // create the atom and short range interaction arrays
232    
233 <  the_atoms = new Atom*[tot_atoms];
234 <  the_molecules = new Molecule[tot_nmol];
233 > #ifdef IS_MPI
234 >
235 >  // divide the molecules among processors here.
236    
237 +  mpiSimulation* mpiSim = new mpiSimulation( simnfo );
238    
239 <  if( tot_SRI ){
240 <    the_sris = new SRI*[tot_SRI];
241 <    the_excludes = new ex_pair[tot_SRI];
239 >  mpiSim->divideLabor();
240 >
241 >  // set up the local variables
242 >  
243 >  int localMol, allMol;
244 >  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
245 >  
246 >  allMol = 0;
247 >  localMol = 0;
248 >  local_atoms = 0;
249 >  local_bonds = 0;
250 >  local_bends = 0;
251 >  local_torsions = 0;
252 >  for( i=0; i<n_components; i++ ){
253 >
254 >    for( j=0; j<components_nmol[i]; j++ ){
255 >      
256 >      if( mpiSim->getMyMolStart() <= allMol &&
257 >          allMol <= mpiSim->getMyMolEnd() ){
258 >        
259 >        local_atoms +=    comp_stamps[i]->getNAtoms();
260 >        local_bonds +=    comp_stamps[i]->getNBonds();
261 >        local_bends +=    comp_stamps[i]->getNBends();
262 >        local_torsions += comp_stamps[i]->getNTorsions();
263 >        localMol++;
264 >      }      
265 >      allMol++;
266 >    }
267    }
268 +  local_SRI = local_bonds + local_bends + local_torsions;
269 +  
270  
271 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
272 +  
273 +  if( local_atoms != simnfo->n_atoms ){
274 +    sprintf( painCave.errMsg,
275 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
276 +             " localAtom (%d) are note equal.\n",
277 +             simnfo->n_atoms,
278 +             local_atoms );
279 +    painCave.isFatal = 1;
280 +    simError();
281 +  }
282 +
283 +  simnfo->n_bonds = local_bonds;
284 +  simnfo->n_bends = local_bends;
285 +  simnfo->n_torsions = local_torsions;
286 +  simnfo->n_SRI = local_SRI;
287 +  simnfo->n_mol = localMol;
288 +
289 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
290 +  MPIcheckPoint();
291 +  
292 +  
293 + #endif // is_mpi
294 +  
295 +
296 +  // create the atom and short range interaction arrays
297 +
298 +  Atom::createArrays(simnfo->n_atoms);
299 +  the_atoms = new Atom*[simnfo->n_atoms];
300 +  the_molecules = new Molecule[simnfo->n_mol];
301 +
302 +
303 +  if( simnfo->n_SRI ){
304 +    the_sris = new SRI*[simnfo->n_SRI];
305 +    the_excludes = new ex_pair[simnfo->n_SRI];
306 +  }
307 +
308    // set the arrays into the SimInfo object
309  
310    simnfo->atoms = the_atoms;
311    simnfo->sr_interactions = the_sris;
312    simnfo->n_exclude = tot_SRI;
313    simnfo->excludes = the_excludes;
145 –  
314  
315 +
316    // initialize the arrays
317 <  
317 >
318    the_ff->setSimInfo( simnfo );
319 <    
319 >
320    makeAtoms();
321  
322    if( tot_bonds ){
# Line 162 | Line 331 | void SimSetup::createSim( void ){
331      makeTorsions();
332    }
333  
165 –  //  makeMolecules();
334  
335    // get some of the tricky things that may still be in the globals
336  
337    if( simnfo->n_dipoles ){
338  
339      if( !the_globals->haveRRF() ){
340 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
341 <      exit(8);
340 >      sprintf( painCave.errMsg,
341 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
342 >      painCave.isFatal = 1;
343 >      simError();
344      }
345      if( !the_globals->haveDielectric() ){
346 <      std::cerr << "SimSetup Error, system has dipoles, but no"
347 <                << " dielectric was set.\n";
348 <      exit(8);
346 >      sprintf( painCave.errMsg,
347 >               "SimSetup Error, system has dipoles, but no"
348 >               " dielectric was set.\n" );
349 >      painCave.isFatal = 1;
350 >      simError();
351      }
352  
353      simnfo->rRF        = the_globals->getRRF();
354      simnfo->dielectric = the_globals->getDielectric();
355    }
356  
357 + #ifdef IS_MPI
358 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
359 +  MPIcheckPoint();
360 + #endif // is_mpi
361 +  
362    if( the_globals->haveBox() ){
363      simnfo->box_x = the_globals->getBox();
364      simnfo->box_y = the_globals->getBox();
365      simnfo->box_z = the_globals->getBox();
366    }
367    else if( the_globals->haveDensity() ){
368 <    
368 >
369      double vol;
370      vol = (double)tot_nmol / the_globals->getDensity();
371      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 197 | Line 374 | void SimSetup::createSim( void ){
374    }
375    else{
376      if( !the_globals->haveBoxX() ){
377 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
378 <      exit(8);
377 >      sprintf( painCave.errMsg,
378 >               "SimSetup error, no periodic BoxX size given.\n" );
379 >      painCave.isFatal = 1;
380 >      simError();
381      }
382      simnfo->box_x = the_globals->getBoxX();
383  
384      if( !the_globals->haveBoxY() ){
385 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
386 <      exit(8);
385 >      sprintf( painCave.errMsg,
386 >               "SimSetup error, no periodic BoxY size given.\n" );
387 >      painCave.isFatal = 1;
388 >      simError();
389      }
390      simnfo->box_y = the_globals->getBoxY();
391  
392      if( !the_globals->haveBoxZ() ){
393 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
394 <      exit(8);
393 >      sprintf( painCave.errMsg,
394 >               "SimSetup error, no periodic BoxZ size given.\n" );
395 >      painCave.isFatal = 1;
396 >      simError();
397      }
398      simnfo->box_z = the_globals->getBoxZ();
399    }
217 –    
218 –  if( the_globals->haveInitialConfig() ){
219 –    InitializeFromFile* fileInit;
220 –    fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
221 –    
222 –    fileInit->read_xyz( simnfo ); // default velocities on
400  
401 <    delete fileInit;    
402 <  }
403 <  else{
404 <    initFromBass();
228 <  }
401 > #ifdef IS_MPI
402 >  strcpy( checkPointMsg, "Box size set up" );
403 >  MPIcheckPoint();
404 > #endif // is_mpi
405  
406 <  if( the_globals->haveFinalConfig() ){
407 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
408 <  }
409 <  else{
410 <    strcpy( simnfo->finalName, inFileName );
411 <    char* endTest;
412 <    int nameLength = strlen( simnfo->finalName );
413 <    endTest = &(simnfo->finalName[nameLength - 5]);
414 <    if( !strcmp( endTest, ".bass" ) ){
415 <      strcpy( endTest, ".eor" );
406 >
407 >
408 > //   if( the_globals->haveInitialConfig() ){
409 > //        InitializeFromFile* fileInit;
410 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
411 >
412 > //     fileInit->read_xyz( simnfo ); // default velocities on
413 >
414 > //     delete fileInit;
415 > //   }
416 > //   else{
417 >
418 > #ifdef IS_MPI
419 >
420 >  // no init from bass
421 >  
422 >  sprintf( painCave.errMsg,
423 >           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
424 >  painCave.isFatal;
425 >  simError();
426 >  
427 > #else
428 >
429 >  initFromBass();
430 >
431 > #endif // is_mpi
432 >
433 > #ifdef IS_MPI
434 >  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
435 >  MPIcheckPoint();
436 > #endif // is_mpi
437 >
438 >
439 >  
440 >
441 >  
442 >  //   }
443 >  
444 > #ifdef IS_MPI
445 >  if( worldRank == 0 ){
446 > #endif // is_mpi
447 >    
448 >    if( the_globals->haveFinalConfig() ){
449 >      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
450      }
241 –    else if( !strcmp( endTest, ".BASS" ) ){
242 –      strcpy( endTest, ".eor" );
243 –    }
451      else{
452 <      endTest = &(simnfo->finalName[nameLength - 4]);
453 <      if( !strcmp( endTest, ".bss" ) ){
452 >      strcpy( simnfo->finalName, inFileName );
453 >      char* endTest;
454 >      int nameLength = strlen( simnfo->finalName );
455 >      endTest = &(simnfo->finalName[nameLength - 5]);
456 >      if( !strcmp( endTest, ".bass" ) ){
457          strcpy( endTest, ".eor" );
458        }
459 <      else if( !strcmp( endTest, ".mdl" ) ){
459 >      else if( !strcmp( endTest, ".BASS" ) ){
460          strcpy( endTest, ".eor" );
461        }
462        else{
463 <        strcat( simnfo->finalName, ".eor" );
463 >        endTest = &(simnfo->finalName[nameLength - 4]);
464 >        if( !strcmp( endTest, ".bss" ) ){
465 >          strcpy( endTest, ".eor" );
466 >        }
467 >        else if( !strcmp( endTest, ".mdl" ) ){
468 >          strcpy( endTest, ".eor" );
469 >        }
470 >        else{
471 >          strcat( simnfo->finalName, ".eor" );
472 >        }
473        }
474      }
256 –  }
475      
476 <  // make the sample and status out names
477 <
478 <  strcpy( simnfo->sampleName, inFileName );
479 <  char* endTest;
480 <  int nameLength = strlen( simnfo->sampleName );
481 <  endTest = &(simnfo->sampleName[nameLength - 5]);
482 <  if( !strcmp( endTest, ".bass" ) ){
265 <    strcpy( endTest, ".dump" );
266 <  }
267 <  else if( !strcmp( endTest, ".BASS" ) ){
268 <    strcpy( endTest, ".dump" );
269 <  }
270 <  else{
271 <    endTest = &(simnfo->sampleName[nameLength - 4]);
272 <    if( !strcmp( endTest, ".bss" ) ){
476 >    // make the sample and status out names
477 >    
478 >    strcpy( simnfo->sampleName, inFileName );
479 >    char* endTest;
480 >    int nameLength = strlen( simnfo->sampleName );
481 >    endTest = &(simnfo->sampleName[nameLength - 5]);
482 >    if( !strcmp( endTest, ".bass" ) ){
483        strcpy( endTest, ".dump" );
484      }
485 <    else if( !strcmp( endTest, ".mdl" ) ){
485 >    else if( !strcmp( endTest, ".BASS" ) ){
486        strcpy( endTest, ".dump" );
487      }
488      else{
489 <      strcat( simnfo->sampleName, ".dump" );
489 >      endTest = &(simnfo->sampleName[nameLength - 4]);
490 >      if( !strcmp( endTest, ".bss" ) ){
491 >        strcpy( endTest, ".dump" );
492 >      }
493 >      else if( !strcmp( endTest, ".mdl" ) ){
494 >        strcpy( endTest, ".dump" );
495 >      }
496 >      else{
497 >        strcat( simnfo->sampleName, ".dump" );
498 >      }
499      }
500 <  }
501 <  
502 <  strcpy( simnfo->statusName, inFileName );
503 <  nameLength = strlen( simnfo->statusName );
504 <  endTest = &(simnfo->statusName[nameLength - 5]);
286 <  if( !strcmp( endTest, ".bass" ) ){
287 <    strcpy( endTest, ".stat" );
288 <  }
289 <  else if( !strcmp( endTest, ".BASS" ) ){
290 <    strcpy( endTest, ".stat" );
291 <  }
292 <  else{
293 <    endTest = &(simnfo->statusName[nameLength - 4]);
294 <    if( !strcmp( endTest, ".bss" ) ){
500 >    
501 >    strcpy( simnfo->statusName, inFileName );
502 >    nameLength = strlen( simnfo->statusName );
503 >    endTest = &(simnfo->statusName[nameLength - 5]);
504 >    if( !strcmp( endTest, ".bass" ) ){
505        strcpy( endTest, ".stat" );
506      }
507 <    else if( !strcmp( endTest, ".mdl" ) ){
507 >    else if( !strcmp( endTest, ".BASS" ) ){
508        strcpy( endTest, ".stat" );
509      }
510      else{
511 <      strcat( simnfo->statusName, ".stat" );
511 >      endTest = &(simnfo->statusName[nameLength - 4]);
512 >      if( !strcmp( endTest, ".bss" ) ){
513 >        strcpy( endTest, ".stat" );
514 >      }
515 >      else if( !strcmp( endTest, ".mdl" ) ){
516 >        strcpy( endTest, ".stat" );
517 >      }
518 >      else{
519 >        strcat( simnfo->statusName, ".stat" );
520 >      }
521      }
522 +    
523 + #ifdef IS_MPI
524    }
525 + #endif // is_mpi
526    
527    // set the status, sample, and themal kick times
528 <
528 >  
529    if( the_globals->haveSampleTime() ){
530 <    simnfo->sampleTime = the_globals->getSampleTime();
530 >    simnfo->sampleTime = the_globals->getSampleTime();
531      simnfo->statusTime = simnfo->sampleTime;
532      simnfo->thermalTime = simnfo->sampleTime;
533    }
534    else{
535 <    simnfo->sampleTime = the_globals->getRunTime();
535 >    simnfo->sampleTime = the_globals->getRunTime();
536      simnfo->statusTime = simnfo->sampleTime;
537      simnfo->thermalTime = simnfo->sampleTime;
538    }
# Line 326 | Line 548 | void SimSetup::createSim( void ){
548    // check for the temperature set flag
549  
550    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
551 <  
552 <    
551 >
552 >
553    // make the longe range forces and the integrator
554 <  
554 >
555    new AllLong( simnfo );
556 <      
556 >
557    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
558    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
559    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
560   }
561  
562   void SimSetup::makeAtoms( void ){
563 <  
563 >
564    int i, j, k, index;
565    double ux, uy, uz, uSqr, u;
566    AtomStamp* current_atom;
567    DirectionalAtom* dAtom;
568 +  int molIndex, molStart, molEnd, nMemb, lMolIndex;
569  
570 +  lMolIndex = 0;
571 +  molIndex = 0;
572    index = 0;
573    for( i=0; i<n_components; i++ ){
574 <    
574 >
575      for( j=0; j<components_nmol[i]; j++ ){
351 –      
352 –      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
353 –        
354 –        current_atom = comp_stamps[i]->getAtom( k );
355 –        if( current_atom->haveOrientation() ){
576  
577 <          dAtom = new DirectionalAtom;
578 <          simnfo->n_oriented++;
579 <          the_atoms[index] = dAtom;
580 <      
581 <          ux = current_atom->getOrntX();
582 <          uy = current_atom->getOrntY();
583 <          uz = current_atom->getOrntZ();
577 > #ifdef IS_MPI
578 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
579 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
580 > #endif // is_mpi        
581 >
582 >        molStart = index;
583 >        nMemb = comp_stamps[i]->getNAtoms();
584 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
585            
586 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
586 >          current_atom = comp_stamps[i]->getAtom( k );
587 >          if( current_atom->haveOrientation() ){
588 >            
589 >            dAtom = new DirectionalAtom(index);
590 >            simnfo->n_oriented++;
591 >            the_atoms[index] = dAtom;
592 >            
593 >            ux = current_atom->getOrntX();
594 >            uy = current_atom->getOrntY();
595 >            uz = current_atom->getOrntZ();
596 >            
597 >            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
598 >            
599 >            u = sqrt( uSqr );
600 >            ux = ux / u;
601 >            uy = uy / u;
602 >            uz = uz / u;
603 >            
604 >            dAtom->setSUx( ux );
605 >            dAtom->setSUy( uy );
606 >            dAtom->setSUz( uz );
607 >          }
608 >          else{
609 >            the_atoms[index] = new GeneralAtom(index);
610 >          }
611 >          the_atoms[index]->setType( current_atom->getType() );
612 >          the_atoms[index]->setIndex( index );
613            
614 <          u = sqrt( uSqr );
615 <          ux = ux / u;
369 <          uy = uy / u;
370 <          uz = uz / u;
371 <          
372 <          dAtom->setSUx( ux );
373 <          dAtom->setSUy( uy );
374 <          dAtom->setSUz( uz );
614 >          // increment the index and repeat;
615 >          index++;
616          }
376 –        else{
377 –          the_atoms[index] = new GeneralAtom;
378 –        }
379 –        the_atoms[index]->setType( current_atom->getType() );
380 –        the_atoms[index]->setIndex( index );
617          
618 <        // increment the index and repeat;
619 <        index++;
618 >        molEnd = index -1;
619 >        the_molecules[lMolIndex].setNMembers( nMemb );
620 >        the_molecules[lMolIndex].setStartAtom( molStart );
621 >        the_molecules[lMolIndex].setEndAtom( molEnd );
622 >        the_molecules[lMolIndex].setStampID( i );
623 >        lMolIndex++;
624 >
625 > #ifdef IS_MPI
626        }
627 + #endif //is_mpi
628 +      
629 +      molIndex++;
630      }
631    }
632 <  
632 >
633    the_ff->initializeAtoms();
634   }
635  
636   void SimSetup::makeBonds( void ){
637  
638 <  int i, j, k, index, offset;
638 >  int i, j, k, index, offset, molIndex;
639    bond_pair* the_bonds;
640    BondStamp* current_bond;
641  
642    the_bonds = new bond_pair[tot_bonds];
643    index = 0;
644    offset = 0;
645 +  molIndex = 0;
646    for( i=0; i<n_components; i++ ){
647 <    
647 >
648      for( j=0; j<components_nmol[i]; j++ ){
403 –      
404 –      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
405 –        
406 –        current_bond = comp_stamps[i]->getBond( k );
407 –        the_bonds[index].a = current_bond->getA() + offset;
408 –        the_bonds[index].b = current_bond->getB() + offset;
649  
650 <        the_excludes[index].i = the_bonds[index].a;
651 <        the_excludes[index].j = the_bonds[index].b;
652 <
653 <        // increment the index and repeat;
654 <        index++;
655 <      }
656 <      offset += comp_stamps[i]->getNAtoms();
657 <    }
650 > #ifdef IS_MPI
651 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
652 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
653 > #endif // is_mpi        
654 >        
655 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
656 >          
657 >          current_bond = comp_stamps[i]->getBond( k );
658 >          the_bonds[index].a = current_bond->getA() + offset;
659 >          the_bonds[index].b = current_bond->getB() + offset;
660 >          
661 >          the_excludes[index].i = the_bonds[index].a;
662 >          the_excludes[index].j = the_bonds[index].b;
663 >          
664 >          // increment the index and repeat;
665 >          index++;
666 >        }
667 >        offset += comp_stamps[i]->getNAtoms();
668 >        
669 > #ifdef IS_MPI
670 >      }
671 > #endif is_mpi
672 >      
673 >      molIndex++;
674 >    }      
675    }
676 <  
676 >
677    the_ff->initializeBonds( the_bonds );
678   }
679  
680   void SimSetup::makeBends( void ){
681  
682 <  int i, j, k, index, offset;
682 >  int i, j, k, index, offset, molIndex;
683    bend_set* the_bends;
684    BendStamp* current_bend;
685  
686    the_bends = new bend_set[tot_bends];
687    index = 0;
688    offset = 0;
689 +  molIndex = 0;
690    for( i=0; i<n_components; i++ ){
691 <    
691 >
692      for( j=0; j<components_nmol[i]; j++ ){
435 –      
436 –      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
437 –        
438 –        current_bend = comp_stamps[i]->getBend( k );
439 –        the_bends[index].a = current_bend->getA() + offset;
440 –        the_bends[index].b = current_bend->getB() + offset;
441 –        the_bends[index].c = current_bend->getC() + offset;
693  
694 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
695 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
694 > #ifdef IS_MPI
695 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
696 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
697 > #endif // is_mpi        
698  
699 <        // increment the index and repeat;
700 <        index++;
699 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
700 >          
701 >          current_bend = comp_stamps[i]->getBend( k );
702 >          the_bends[index].a = current_bend->getA() + offset;
703 >          the_bends[index].b = current_bend->getB() + offset;
704 >          the_bends[index].c = current_bend->getC() + offset;
705 >          
706 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
707 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
708 >          
709 >          // increment the index and repeat;
710 >          index++;
711 >        }
712 >        offset += comp_stamps[i]->getNAtoms();
713 >        
714 > #ifdef IS_MPI
715        }
716 <      offset += comp_stamps[i]->getNAtoms();
716 > #endif //is_mpi
717 >
718 >      molIndex++;
719      }
720    }
721 <  
721 >
722    the_ff->initializeBends( the_bends );
723   }
724  
725   void SimSetup::makeTorsions( void ){
726  
727 <  int i, j, k, index, offset;
727 >  int i, j, k, index, offset, molIndex;
728    torsion_set* the_torsions;
729    TorsionStamp* current_torsion;
730  
731    the_torsions = new torsion_set[tot_torsions];
732    index = 0;
733    offset = 0;
734 +  molIndex = 0;
735    for( i=0; i<n_components; i++ ){
736 <    
736 >
737      for( j=0; j<components_nmol[i]; j++ ){
738 <      
738 >
739 > #ifdef IS_MPI
740 >      if( simnfo->mpiSim->getMyMolStart() <= molIndex &&
741 >          molIndex <= simnfo->mpiSim->getMyMolEnd() ){
742 > #endif // is_mpi        
743 >
744        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
745 <        
745 >
746          current_torsion = comp_stamps[i]->getTorsion( k );
747          the_torsions[index].a = current_torsion->getA() + offset;
748          the_torsions[index].b = current_torsion->getB() + offset;
# Line 481 | Line 756 | void SimSetup::makeTorsions( void ){
756          index++;
757        }
758        offset += comp_stamps[i]->getNAtoms();
759 +
760 + #ifdef IS_MPI
761 +      }
762 + #endif //is_mpi      
763 +
764 +      molIndex++;
765      }
766    }
486 –  
487 –  the_ff->initializeTorsions( the_torsions );
488 – }
767  
768 < void SimSetup::makeMolecules( void ){
491 <
492 <  int i,j,k;
493 <  
494 <  for( i=0; i<n_components; i++ ){
495 <    
496 <    for( j=0; j<components_nmol[i]; j++ ){
497 <      
498 <      
499 <
768 >  the_ff->initializeTorsions( the_torsions );
769   }
770  
771   void SimSetup::initFromBass( void ){
# Line 526 | Line 795 | void SimSetup::initFromBass( void ){
795      n_per_extra = (int)ceil( temp1 );
796  
797      if( n_per_extra > 4){
798 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
799 <      exit(8);
798 >      sprintf( painCave.errMsg,
799 >               "SimSetup error. There has been an error in constructing"
800 >               " the non-complete lattice.\n" );
801 >      painCave.isFatal = 1;
802 >      simError();
803      }
804    }
805    else{
# Line 536 | Line 808 | void SimSetup::initFromBass( void ){
808      celly = simnfo->box_y / temp3;
809      cellz = simnfo->box_z / temp3;
810    }
811 <  
811 >
812    current_mol = 0;
813    current_comp_mol = 0;
814    current_comp = 0;
815    current_atom_ndx = 0;
816 <  
816 >
817    for( i=0; i < n_cells ; i++ ){
818      for( j=0; j < n_cells; j++ ){
819        for( k=0; k < n_cells; k++ ){
820 <        
820 >
821          makeElement( i * cellx,
822                       j * celly,
823                       k * cellz );
824 <        
824 >
825          makeElement( i * cellx + 0.5 * cellx,
826                       j * celly + 0.5 * celly,
827                       k * cellz );
828 <        
828 >
829          makeElement( i * cellx,
830                       j * celly + 0.5 * celly,
831                       k * cellz + 0.5 * cellz );
832 <        
832 >
833          makeElement( i * cellx + 0.5 * cellx,
834                       j * celly,
835                       k * cellz + 0.5 * cellz );
# Line 567 | Line 839 | void SimSetup::initFromBass( void ){
839  
840    if( have_extra ){
841      done = 0;
842 <    
842 >
843      int start_ndx;
844      for( i=0; i < (n_cells+1) && !done; i++ ){
845        for( j=0; j < (n_cells+1) && !done; j++ ){
846 <        
846 >
847          if( i < n_cells ){
848 <          
848 >
849            if( j < n_cells ){
850              start_ndx = n_cells;
851            }
852            else start_ndx = 0;
853          }
854          else start_ndx = 0;
855 <        
855 >
856          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
857 <          
857 >
858            makeElement( i * cellx,
859                         j * celly,
860                         k * cellz );
861            done = ( current_mol >= tot_nmol );
862 <          
862 >
863            if( !done && n_per_extra > 1 ){
864              makeElement( i * cellx + 0.5 * cellx,
865                           j * celly + 0.5 * celly,
866                           k * cellz );
867              done = ( current_mol >= tot_nmol );
868            }
869 <          
869 >
870            if( !done && n_per_extra > 2){
871              makeElement( i * cellx,
872                           j * celly + 0.5 * celly,
873                           k * cellz + 0.5 * cellz );
874              done = ( current_mol >= tot_nmol );
875            }
876 <          
876 >
877            if( !done && n_per_extra > 3){
878              makeElement( i * cellx + 0.5 * cellx,
879                           j * celly,
# Line 612 | Line 884 | void SimSetup::initFromBass( void ){
884        }
885      }
886    }
887 <  
888 <  
887 >
888 >
889    for( i=0; i<simnfo->n_atoms; i++ ){
890      simnfo->atoms[i]->set_vx( 0.0 );
891      simnfo->atoms[i]->set_vy( 0.0 );
# Line 629 | Line 901 | void SimSetup::makeElement( double x, double y, double
901    double rotMat[3][3];
902  
903    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
904 <    
904 >
905      current_atom = comp_stamps[current_comp]->getAtom( k );
906      if( !current_atom->havePosition() ){
907 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
908 <                << ", atom " << current_atom->getType()
909 <                << " does not have a position specified.\n"
910 <                << "The initialization routine is unable to give a start"
911 <                << " position.\n";
912 <      exit(8);
907 >      sprintf( painCave.errMsg,
908 >               "SimSetup:initFromBass error.\n"
909 >               "\tComponent %s, atom %s does not have a position specified.\n"
910 >               "\tThe initialization routine is unable to give a start"
911 >               " position.\n",
912 >               comp_stamps[current_comp]->getID(),
913 >               current_atom->getType() );
914 >      painCave.isFatal = 1;
915 >      simError();
916      }
917 <    
917 >
918      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
919      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
920      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
921 <    
921 >
922      if( the_atoms[current_atom_ndx]->isDirectional() ){
923 <      
923 >
924        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
925 <      
925 >
926        rotMat[0][0] = 1.0;
927        rotMat[0][1] = 0.0;
928        rotMat[0][2] = 0.0;
# Line 665 | Line 940 | void SimSetup::makeElement( double x, double y, double
940  
941      current_atom_ndx++;
942    }
943 <  
943 >
944    current_mol++;
945    current_comp_mol++;
946  
947    if( current_comp_mol >= components_nmol[current_comp] ){
948 <    
948 >
949      current_comp_mol = 0;
950      current_comp++;
951    }

Diff Legend

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