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 176 by mmeineke, Thu Nov 14 22:00:44 2002 UTC vs.
Revision 261 by chuckv, Mon Feb 3 21:15:59 2003 UTC

# Line 10 | Line 10
10  
11   #ifdef IS_MPI
12   #include "mpiBASS.h"
13 + #include "mpiSimulation.hpp"
14   #include "bassDiag.hpp"
15   #endif
16  
# Line 73 | Line 74 | void SimSetup::createSim( void ){
74  
75    MakeStamps *the_stamps;
76    Globals* the_globals;
77 <  int i;
77 >  int i, j;
78  
79    // get the stamps and globals;
80    the_stamps = stamps;
# Line 92 | Line 93 | void SimSetup::createSim( void ){
93    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
94    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
95    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
96 +  else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
97    else{
98      sprintf( painCave.errMsg,
99               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 119 | Line 121 | void SimSetup::createSim( void ){
121  
122        if( !the_components[i]->haveNMol() ){
123          // we have a problem
124 <        std::cerr << "SimSetup Error. No global NMol or component NMol"
125 <                  << " given. Cannot calculate the number of atoms.\n";
126 <        exit( 8 );
124 >        sprintf( painCave.errMsg,
125 >                 "SimSetup Error. No global NMol or component NMol"
126 >                 " given. Cannot calculate the number of atoms.\n" );
127 >        painCave.isFatal = 1;
128 >        simError();
129        }
130  
131        tot_nmol += the_components[i]->getNMol();
# Line 129 | Line 133 | void SimSetup::createSim( void ){
133      }
134    }
135    else{
136 <    std::cerr << "NOT A SUPPORTED FEATURE\n";
137 <
138 < //     tot_nmol = the_globals->getNMol();
139 <
140 < //     //we have the total number of molecules, now we check for molfractions
141 < //     for( i=0; i<n_components; i++ ){
142 <
143 < //       if( !the_components[i]->haveMolFraction() ){
144 <
145 < //      if( !the_components[i]->haveNMol() ){
146 < //        //we have a problem
147 < //        std::cerr << "SimSetup error. Neither molFraction nor "
148 < //                  << " nMol was given in component
149 <
136 >    sprintf( painCave.errMsg,
137 >             "SimSetup error.\n"
138 >             "\tSorry, the ability to specify total"
139 >             " nMols and then give molfractions in the components\n"
140 >             "\tis not currently supported."
141 >             " Please give nMol in the components.\n" );
142 >    painCave.isFatal = 1;
143 >    simError();
144 >    
145 >    
146 >    //     tot_nmol = the_globals->getNMol();
147 >    
148 >    //   //we have the total number of molecules, now we check for molfractions
149 >    //     for( i=0; i<n_components; i++ ){
150 >    
151 >    //       if( !the_components[i]->haveMolFraction() ){
152 >    
153 >    //  if( !the_components[i]->haveNMol() ){
154 >    //    //we have a problem
155 >    //    std::cerr << "SimSetup error. Neither molFraction nor "
156 >    //              << " nMol was given in component
157 >    
158    }
159  
160 + #ifdef IS_MPI
161 +  strcpy( checkPointMsg, "Have the number of components" );
162 +  MPIcheckPoint();
163 + #endif // is_mpi
164 +
165    // make an array of molecule stamps that match the components used.
166 +  // also extract the used stamps out into a separate linked list
167  
168 +  simnfo->nComponents = n_components;
169 +  simnfo->componentsNmol = components_nmol;
170 +  simnfo->compStamps = comp_stamps;
171 +  simnfo->headStamp = new LinkedMolStamp();
172 +  
173 +  char* id;
174 +  LinkedMolStamp* headStamp = simnfo->headStamp;
175 +  LinkedMolStamp* currentStamp = NULL;
176    for( i=0; i<n_components; i++ ){
177  
178 <    comp_stamps[i] =
179 <      the_stamps->getMolecule( the_components[i]->getType() );
178 >    id = the_components[i]->getType();
179 >    comp_stamps[i] = NULL;
180 >    
181 >    // check to make sure the component isn't already in the list
182 >
183 >    comp_stamps[i] = headStamp->match( id );
184 >    if( comp_stamps[i] == NULL ){
185 >      
186 >      // extract the component from the list;
187 >      
188 >      currentStamp = the_stamps->extractMolStamp( id );
189 >      if( currentStamp == NULL ){
190 >        sprintf( painCave.errMsg,
191 >                 "SimSetup error: Component \"%s\" was not found in the "
192 >                 "list of declared molecules\n",
193 >                 id );
194 >        painCave.isFatal = 1;
195 >        simError();
196 >      }
197 >      
198 >      headStamp->add( currentStamp );
199 >      comp_stamps[i] = headStamp->match( id );
200 >    }
201    }
202  
203 + #ifdef IS_MPI
204 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
205 +  MPIcheckPoint();
206 + #endif // is_mpi
207 +  
208  
209  
210 +
211    // caclulate the number of atoms, bonds, bends and torsions
212  
213    tot_atoms = 0;
# Line 162 | Line 215 | void SimSetup::createSim( void ){
215    tot_bends = 0;
216    tot_torsions = 0;
217    for( i=0; i<n_components; i++ ){
218 <
219 <    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
220 <    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
221 <    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
218 >    
219 >    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
220 >    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
221 >    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
222      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
223    }
224  
# Line 176 | Line 229 | void SimSetup::createSim( void ){
229    simnfo->n_bends = tot_bends;
230    simnfo->n_torsions = tot_torsions;
231    simnfo->n_SRI = tot_SRI;
232 +  simnfo->n_mol = tot_nmol;
233 +
234 +  
235 + #ifdef IS_MPI
236  
237 <  // create the atom and short range interaction arrays
237 >  // divide the molecules among processors here.
238 >  
239 >  mpiSim = new mpiSimulation( simnfo );
240 >  
241 >  
242  
243 <  the_atoms = new Atom*[tot_atoms];
183 <  Atom::createArrays(tot_atoms);
184 <  the_molecules = new Molecule[tot_nmol];
243 >  globalIndex = mpiSim->divideLabor();
244  
245  
187 –  if( tot_SRI ){
188 –    the_sris = new SRI*[tot_SRI];
189 –    the_excludes = new ex_pair[tot_SRI];
190 –  }
246  
247 <  // set the arrays into the SimInfo object
247 >  // set up the local variables
248 >  
249 >  int localMol, allMol;
250 >  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
251 >  
252 >  allMol = 0;
253 >  localMol = 0;
254 >  local_atoms = 0;
255 >  local_bonds = 0;
256 >  local_bends = 0;
257 >  local_torsions = 0;
258 >  for( i=0; i<n_components; i++ ){
259  
260 <  simnfo->atoms = the_atoms;
261 <  simnfo->sr_interactions = the_sris;
262 <  simnfo->n_exclude = tot_SRI;
263 <  simnfo->excludes = the_excludes;
260 >    for( j=0; j<components_nmol[i]; j++ ){
261 >      
262 >      if( mpiSim->getMyMolStart() <= allMol &&
263 >          allMol <= mpiSim->getMyMolEnd() ){
264 >        
265 >        local_atoms +=    comp_stamps[i]->getNAtoms();
266 >        local_bonds +=    comp_stamps[i]->getNBonds();
267 >        local_bends +=    comp_stamps[i]->getNBends();
268 >        local_torsions += comp_stamps[i]->getNTorsions();
269 >        localMol++;
270 >      }      
271 >      allMol++;
272 >    }
273 >  }
274 >  local_SRI = local_bonds + local_bends + local_torsions;
275 >  
276  
277 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
278 +  
279 +  if( local_atoms != simnfo->n_atoms ){
280 +    sprintf( painCave.errMsg,
281 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
282 +             " localAtom (%d) are note equal.\n",
283 +             simnfo->n_atoms,
284 +             local_atoms );
285 +    painCave.isFatal = 1;
286 +    simError();
287 +  }
288  
289 <  // initialize the arrays
289 >  simnfo->n_bonds = local_bonds;
290 >  simnfo->n_bends = local_bends;
291 >  simnfo->n_torsions = local_torsions;
292 >  simnfo->n_SRI = local_SRI;
293 >  simnfo->n_mol = localMol;
294  
295 <  the_ff->setSimInfo( simnfo );
295 >  strcpy( checkPointMsg, "Passed nlocal consistency check." );
296 >  MPIcheckPoint();
297 >  
298 >  
299 > #endif // is_mpi
300 >  
301  
302 <  makeAtoms();
302 >  // create the atom and short range interaction arrays
303  
304 <  if( tot_bonds ){
305 <    makeBonds();
306 <  }
304 >  Atom::createArrays(simnfo->n_atoms);
305 >  the_atoms = new Atom*[simnfo->n_atoms];
306 >  the_molecules = new Molecule[simnfo->n_mol];
307  
308 <  if( tot_bends ){
309 <    makeBends();
310 <  }
311 <
214 <  if( tot_torsions ){
215 <    makeTorsions();
308 >
309 >  if( simnfo->n_SRI ){
310 >    the_sris = new SRI*[simnfo->n_SRI];
311 >    the_excludes = new ex_pair[simnfo->n_SRI];
312    }
313  
314 <  //  makeMolecules();
314 >  // set the arrays into the SimInfo object
315  
316 +  simnfo->atoms = the_atoms;
317 +  simnfo->sr_interactions = the_sris;
318 +  simnfo->n_exclude = tot_SRI;
319 +  simnfo->excludes = the_excludes;
320 +
321 +
322    // get some of the tricky things that may still be in the globals
323  
324    if( simnfo->n_dipoles ){
325  
326      if( !the_globals->haveRRF() ){
327 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
328 <      exit(8);
327 >      sprintf( painCave.errMsg,
328 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
329 >      painCave.isFatal = 1;
330 >      simError();
331      }
332      if( !the_globals->haveDielectric() ){
333 <      std::cerr << "SimSetup Error, system has dipoles, but no"
334 <                << " dielectric was set.\n";
335 <      exit(8);
333 >      sprintf( painCave.errMsg,
334 >               "SimSetup Error, system has dipoles, but no"
335 >               " dielectric was set.\n" );
336 >      painCave.isFatal = 1;
337 >      simError();
338      }
339  
340      simnfo->rRF        = the_globals->getRRF();
341      simnfo->dielectric = the_globals->getDielectric();
342    }
343  
344 + #ifdef IS_MPI
345 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
346 +  MPIcheckPoint();
347 + #endif // is_mpi
348 +  
349    if( the_globals->haveBox() ){
350      simnfo->box_x = the_globals->getBox();
351      simnfo->box_y = the_globals->getBox();
# Line 250 | Line 361 | void SimSetup::createSim( void ){
361    }
362    else{
363      if( !the_globals->haveBoxX() ){
364 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
365 <      exit(8);
364 >      sprintf( painCave.errMsg,
365 >               "SimSetup error, no periodic BoxX size given.\n" );
366 >      painCave.isFatal = 1;
367 >      simError();
368      }
369      simnfo->box_x = the_globals->getBoxX();
370  
371      if( !the_globals->haveBoxY() ){
372 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
373 <      exit(8);
372 >      sprintf( painCave.errMsg,
373 >               "SimSetup error, no periodic BoxY size given.\n" );
374 >      painCave.isFatal = 1;
375 >      simError();
376      }
377      simnfo->box_y = the_globals->getBoxY();
378  
379      if( !the_globals->haveBoxZ() ){
380 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
381 <      exit(8);
380 >      sprintf( painCave.errMsg,
381 >               "SimSetup error, no periodic BoxZ size given.\n" );
382 >      painCave.isFatal = 1;
383 >      simError();
384      }
385      simnfo->box_z = the_globals->getBoxZ();
386    }
387  
388 + #ifdef IS_MPI
389 +  strcpy( checkPointMsg, "Box size set up" );
390 +  MPIcheckPoint();
391 + #endif // is_mpi
392  
272 – //   if( the_globals->haveInitialConfig() ){
273 – //        InitializeFromFile* fileInit;
274 – //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
393  
394 < //     fileInit->read_xyz( simnfo ); // default velocities on
394 >  // initialize the arrays
395  
396 < //     delete fileInit;
279 < //   }
280 < //   else{
396 >  the_ff->setSimInfo( simnfo );
397  
398 <    initFromBass();
398 >  makeAtoms();
399 >  //
400 >  if( tot_bonds ){
401 >    makeBonds();
402 >  }
403  
404 +  if( tot_bends ){
405 +    makeBends();
406 +  }
407  
408 < //   }
408 >  if( tot_torsions ){
409 >    makeTorsions();
410 >  }
411  
287 – //   if( the_globals->haveFinalConfig() ){
288 – //     strcpy( simnfo->finalName, the_globals->getFinalConfig() );
289 – //   }
290 – //   else{
291 – //     strcpy( simnfo->finalName, inFileName );
292 – //     char* endTest;
293 – //     int nameLength = strlen( simnfo->finalName );
294 – //     endTest = &(simnfo->finalName[nameLength - 5]);
295 – //     if( !strcmp( endTest, ".bass" ) ){
296 – //       strcpy( endTest, ".eor" );
297 – //     }
298 – //     else if( !strcmp( endTest, ".BASS" ) ){
299 – //       strcpy( endTest, ".eor" );
300 – //     }
301 – //     else{
302 – //       endTest = &(simnfo->finalName[nameLength - 4]);
303 – //       if( !strcmp( endTest, ".bss" ) ){
304 – //      strcpy( endTest, ".eor" );
305 – //       }
306 – //       else if( !strcmp( endTest, ".mdl" ) ){
307 – //      strcpy( endTest, ".eor" );
308 – //       }
309 – //       else{
310 – //      strcat( simnfo->finalName, ".eor" );
311 – //       }
312 – //     }
313 – //   }
412  
315 – //   // make the sample and status out names
413  
317 – //   strcpy( simnfo->sampleName, inFileName );
318 – //   char* endTest;
319 – //   int nameLength = strlen( simnfo->sampleName );
320 – //   endTest = &(simnfo->sampleName[nameLength - 5]);
321 – //   if( !strcmp( endTest, ".bass" ) ){
322 – //     strcpy( endTest, ".dump" );
323 – //   }
324 – //   else if( !strcmp( endTest, ".BASS" ) ){
325 – //     strcpy( endTest, ".dump" );
326 – //   }
327 – //   else{
328 – //     endTest = &(simnfo->sampleName[nameLength - 4]);
329 – //     if( !strcmp( endTest, ".bss" ) ){
330 – //       strcpy( endTest, ".dump" );
331 – //     }
332 – //     else if( !strcmp( endTest, ".mdl" ) ){
333 – //       strcpy( endTest, ".dump" );
334 – //     }
335 – //     else{
336 – //       strcat( simnfo->sampleName, ".dump" );
337 – //     }
338 – //   }
414  
340 – //   strcpy( simnfo->statusName, inFileName );
341 – //   nameLength = strlen( simnfo->statusName );
342 – //   endTest = &(simnfo->statusName[nameLength - 5]);
343 – //   if( !strcmp( endTest, ".bass" ) ){
344 – //     strcpy( endTest, ".stat" );
345 – //   }
346 – //   else if( !strcmp( endTest, ".BASS" ) ){
347 – //     strcpy( endTest, ".stat" );
348 – //   }
349 – //   else{
350 – //     endTest = &(simnfo->statusName[nameLength - 4]);
351 – //     if( !strcmp( endTest, ".bss" ) ){
352 – //       strcpy( endTest, ".stat" );
353 – //     }
354 – //     else if( !strcmp( endTest, ".mdl" ) ){
355 – //       strcpy( endTest, ".stat" );
356 – //     }
357 – //     else{
358 – //       strcat( simnfo->statusName, ".stat" );
359 – //     }
360 – //   }
415  
416  
417 <  // set the status, sample, and themal kick times
417 > if( the_globals->haveInitialConfig() ){
418 >
419 >     InitializeFromFile* fileInit;
420 > #ifdef IS_MPI // is_mpi
421 >     if( worldRank == 0 ){
422 > #endif //is_mpi
423 >   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
424 > #ifdef IS_MPI
425 >     }else fileInit = new InitializeFromFile( NULL );
426 > #endif
427 >   fileInit->read_xyz( simnfo ); // default velocities on
428  
429 +   delete fileInit;
430 + }
431 + else{
432 +
433 + #ifdef IS_MPI
434 +
435 +  // no init from bass
436 +  
437 +  sprintf( painCave.errMsg,
438 +           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
439 +  painCave.isFatal;
440 +  simError();
441 +  
442 + #else
443 +
444 +  initFromBass();
445 +
446 +
447 + #endif
448 + }
449 +
450 + #ifdef IS_MPI
451 +  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
452 +  MPIcheckPoint();
453 + #endif // is_mpi
454 +
455 +
456 +  
457 +
458 +  
459 +
460 +  
461 + #ifdef IS_MPI
462 +  if( worldRank == 0 ){
463 + #endif // is_mpi
464 +    
465 +    if( the_globals->haveFinalConfig() ){
466 +      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
467 +    }
468 +    else{
469 +      strcpy( simnfo->finalName, inFileName );
470 +      char* endTest;
471 +      int nameLength = strlen( simnfo->finalName );
472 +      endTest = &(simnfo->finalName[nameLength - 5]);
473 +      if( !strcmp( endTest, ".bass" ) ){
474 +        strcpy( endTest, ".eor" );
475 +      }
476 +      else if( !strcmp( endTest, ".BASS" ) ){
477 +        strcpy( endTest, ".eor" );
478 +      }
479 +      else{
480 +        endTest = &(simnfo->finalName[nameLength - 4]);
481 +        if( !strcmp( endTest, ".bss" ) ){
482 +          strcpy( endTest, ".eor" );
483 +        }
484 +        else if( !strcmp( endTest, ".mdl" ) ){
485 +          strcpy( endTest, ".eor" );
486 +        }
487 +        else{
488 +          strcat( simnfo->finalName, ".eor" );
489 +        }
490 +      }
491 +    }
492 +    
493 +    // make the sample and status out names
494 +    
495 +    strcpy( simnfo->sampleName, inFileName );
496 +    char* endTest;
497 +    int nameLength = strlen( simnfo->sampleName );
498 +    endTest = &(simnfo->sampleName[nameLength - 5]);
499 +    if( !strcmp( endTest, ".bass" ) ){
500 +      strcpy( endTest, ".dump" );
501 +    }
502 +    else if( !strcmp( endTest, ".BASS" ) ){
503 +      strcpy( endTest, ".dump" );
504 +    }
505 +    else{
506 +      endTest = &(simnfo->sampleName[nameLength - 4]);
507 +      if( !strcmp( endTest, ".bss" ) ){
508 +        strcpy( endTest, ".dump" );
509 +      }
510 +      else if( !strcmp( endTest, ".mdl" ) ){
511 +        strcpy( endTest, ".dump" );
512 +      }
513 +      else{
514 +        strcat( simnfo->sampleName, ".dump" );
515 +      }
516 +    }
517 +    
518 +    strcpy( simnfo->statusName, inFileName );
519 +    nameLength = strlen( simnfo->statusName );
520 +    endTest = &(simnfo->statusName[nameLength - 5]);
521 +    if( !strcmp( endTest, ".bass" ) ){
522 +      strcpy( endTest, ".stat" );
523 +    }
524 +    else if( !strcmp( endTest, ".BASS" ) ){
525 +      strcpy( endTest, ".stat" );
526 +    }
527 +    else{
528 +      endTest = &(simnfo->statusName[nameLength - 4]);
529 +      if( !strcmp( endTest, ".bss" ) ){
530 +        strcpy( endTest, ".stat" );
531 +      }
532 +      else if( !strcmp( endTest, ".mdl" ) ){
533 +        strcpy( endTest, ".stat" );
534 +      }
535 +      else{
536 +        strcat( simnfo->statusName, ".stat" );
537 +      }
538 +    }
539 +    
540 + #ifdef IS_MPI
541 +  }
542 + #endif // is_mpi
543 +  
544 +  // set the status, sample, and themal kick times
545 +  
546    if( the_globals->haveSampleTime() ){
547      simnfo->sampleTime = the_globals->getSampleTime();
548      simnfo->statusTime = simnfo->sampleTime;
# Line 386 | Line 567 | void SimSetup::createSim( void ){
567    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
568  
569  
570 <  // make the longe range forces and the integrator
570 > //   // make the longe range forces and the integrator
571  
572 <  new AllLong( simnfo );
572 > //   new AllLong( simnfo );
573  
574 <  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
574 >  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
575    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
576    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
577 +  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
578 +
579   }
580  
581   void SimSetup::makeAtoms( void ){
# Line 401 | Line 584 | void SimSetup::makeAtoms( void ){
584    double ux, uy, uz, uSqr, u;
585    AtomStamp* current_atom;
586    DirectionalAtom* dAtom;
587 <  int molIndex, molStart, molEnd, nMemb;
587 >  int molIndex, molStart, molEnd, nMemb, lMolIndex;
588  
589 <
589 >  lMolIndex = 0;
590    molIndex = 0;
591    index = 0;
592    for( i=0; i<n_components; i++ ){
593  
594      for( j=0; j<components_nmol[i]; j++ ){
595  
596 <      molStart = index;
597 <      nMemb = comp_stamps[i]->getNAtoms();
598 <      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
596 > #ifdef IS_MPI
597 >      if( mpiSim->getMyMolStart() <= molIndex &&
598 >          molIndex <= mpiSim->getMyMolEnd() ){
599 > #endif // is_mpi        
600  
601 <        current_atom = comp_stamps[i]->getAtom( k );
602 <        if( current_atom->haveOrientation() ){
603 <
604 <          dAtom = new DirectionalAtom(index);
605 <          simnfo->n_oriented++;
606 <          the_atoms[index] = dAtom;
607 <
608 <          ux = current_atom->getOrntX();
609 <          uy = current_atom->getOrntY();
610 <          uz = current_atom->getOrntZ();
611 <
612 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
613 <
614 <          u = sqrt( uSqr );
615 <          ux = ux / u;
616 <          uy = uy / u;
617 <          uz = uz / u;
618 <
619 <          dAtom->setSUx( ux );
620 <          dAtom->setSUy( uy );
621 <          dAtom->setSUz( uz );
622 <        }
623 <        else{
624 <          the_atoms[index] = new GeneralAtom(index);
625 <        }
626 <        the_atoms[index]->setType( current_atom->getType() );
627 <        the_atoms[index]->setIndex( index );
628 <
629 <        // increment the index and repeat;
630 <        index++;
631 <      }
601 >        molStart = index;
602 >        nMemb = comp_stamps[i]->getNAtoms();
603 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
604 >          
605 >          current_atom = comp_stamps[i]->getAtom( k );
606 >          if( current_atom->haveOrientation() ){
607 >            
608 >            dAtom = new DirectionalAtom(index);
609 >            simnfo->n_oriented++;
610 >            the_atoms[index] = dAtom;
611 >            
612 >            ux = current_atom->getOrntX();
613 >            uy = current_atom->getOrntY();
614 >            uz = current_atom->getOrntZ();
615 >            
616 >            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
617 >            
618 >            u = sqrt( uSqr );
619 >            ux = ux / u;
620 >            uy = uy / u;
621 >            uz = uz / u;
622 >            
623 >            dAtom->setSUx( ux );
624 >            dAtom->setSUy( uy );
625 >            dAtom->setSUz( uz );
626 >          }
627 >          else{
628 >            the_atoms[index] = new GeneralAtom(index);
629 >          }
630 >          the_atoms[index]->setType( current_atom->getType() );
631 >          the_atoms[index]->setIndex( index );
632 >          
633 >          // increment the index and repeat;
634 >          index++;
635 >        }
636 >        
637 >        molEnd = index -1;
638 >        the_molecules[lMolIndex].setNMembers( nMemb );
639 >        the_molecules[lMolIndex].setStartAtom( molStart );
640 >        the_molecules[lMolIndex].setEndAtom( molEnd );
641 >        the_molecules[lMolIndex].setStampID( i );
642 >        lMolIndex++;
643  
644 <      molEnd = index -1;
645 <      the_molecules[molIndex].setNMembers( nMemb );
646 <      the_molecules[molIndex].setStartAtom( molStart );
647 <      the_molecules[molIndex].setEndAtom( molEnd );
644 > #ifdef IS_MPI
645 >      }
646 > #endif //is_mpi
647 >      
648        molIndex++;
454 –
649      }
650    }
651  
652 + #ifdef IS_MPI
653 +    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
654 +    
655 +    delete[] globalIndex;
656 +
657 +    mpiSim->mpiRefresh();
658 + #endif //IS_MPI
659 +          
660    the_ff->initializeAtoms();
661   }
662  
663   void SimSetup::makeBonds( void ){
664  
665 <  int i, j, k, index, offset;
665 >  int i, j, k, index, offset, molIndex;
666    bond_pair* the_bonds;
667    BondStamp* current_bond;
668  
669    the_bonds = new bond_pair[tot_bonds];
670    index = 0;
671    offset = 0;
672 +  molIndex = 0;
673 +
674    for( i=0; i<n_components; i++ ){
675  
676      for( j=0; j<components_nmol[i]; j++ ){
677  
678 <      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
679 <
680 <        current_bond = comp_stamps[i]->getBond( k );
681 <        the_bonds[index].a = current_bond->getA() + offset;
682 <        the_bonds[index].b = current_bond->getB() + offset;
683 <
684 <        the_excludes[index].i = the_bonds[index].a;
685 <        the_excludes[index].j = the_bonds[index].b;
686 <
687 <        // increment the index and repeat;
688 <        index++;
678 > #ifdef IS_MPI
679 >      if( mpiSim->getMyMolStart() <= molIndex &&
680 >          molIndex <= mpiSim->getMyMolEnd() ){
681 > #endif // is_mpi        
682 >        
683 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
684 >          
685 >          current_bond = comp_stamps[i]->getBond( k );
686 >          the_bonds[index].a = current_bond->getA() + offset;
687 >          the_bonds[index].b = current_bond->getB() + offset;
688 >          
689 >          the_excludes[index].i = the_bonds[index].a;
690 >          the_excludes[index].j = the_bonds[index].b;
691 >          
692 >          // increment the index and repeat;
693 >          index++;
694 >        }
695 >        offset += comp_stamps[i]->getNAtoms();
696 >        
697 > #ifdef IS_MPI
698        }
699 <      offset += comp_stamps[i]->getNAtoms();
700 <    }
699 > #endif //is_mpi
700 >      
701 >      molIndex++;
702 >    }      
703    }
704  
705    the_ff->initializeBonds( the_bonds );
# Line 492 | Line 707 | void SimSetup::makeBends( void ){
707  
708   void SimSetup::makeBends( void ){
709  
710 <  int i, j, k, index, offset;
710 >  int i, j, k, index, offset, molIndex;
711    bend_set* the_bends;
712    BendStamp* current_bend;
713  
714    the_bends = new bend_set[tot_bends];
715    index = 0;
716    offset = 0;
717 +  molIndex = 0;
718    for( i=0; i<n_components; i++ ){
719  
720      for( j=0; j<components_nmol[i]; j++ ){
721  
722 <      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
722 > #ifdef IS_MPI
723 >      if( mpiSim->getMyMolStart() <= molIndex &&
724 >          molIndex <= mpiSim->getMyMolEnd() ){
725 > #endif // is_mpi        
726  
727 <        current_bend = comp_stamps[i]->getBend( k );
728 <        the_bends[index].a = current_bend->getA() + offset;
729 <        the_bends[index].b = current_bend->getB() + offset;
730 <        the_bends[index].c = current_bend->getC() + offset;
731 <
732 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
733 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
734 <
735 <        // increment the index and repeat;
736 <        index++;
727 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
728 >          
729 >          current_bend = comp_stamps[i]->getBend( k );
730 >          the_bends[index].a = current_bend->getA() + offset;
731 >          the_bends[index].b = current_bend->getB() + offset;
732 >          the_bends[index].c = current_bend->getC() + offset;
733 >          
734 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
735 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
736 >          
737 >          // increment the index and repeat;
738 >          index++;
739 >        }
740 >        offset += comp_stamps[i]->getNAtoms();
741 >        
742 > #ifdef IS_MPI
743        }
744 <      offset += comp_stamps[i]->getNAtoms();
744 > #endif //is_mpi
745 >
746 >      molIndex++;
747      }
748    }
749  
# Line 525 | Line 752 | void SimSetup::makeTorsions( void ){
752  
753   void SimSetup::makeTorsions( void ){
754  
755 <  int i, j, k, index, offset;
755 >  int i, j, k, index, offset, molIndex;
756    torsion_set* the_torsions;
757    TorsionStamp* current_torsion;
758  
759    the_torsions = new torsion_set[tot_torsions];
760    index = 0;
761    offset = 0;
762 +  molIndex = 0;
763    for( i=0; i<n_components; i++ ){
764  
765      for( j=0; j<components_nmol[i]; j++ ){
766  
767 + #ifdef IS_MPI
768 +      if( mpiSim->getMyMolStart() <= molIndex &&
769 +          molIndex <= mpiSim->getMyMolEnd() ){
770 + #endif // is_mpi        
771 +
772        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
773  
774          current_torsion = comp_stamps[i]->getTorsion( k );
# Line 551 | Line 784 | void SimSetup::makeTorsions( void ){
784          index++;
785        }
786        offset += comp_stamps[i]->getNAtoms();
787 +
788 + #ifdef IS_MPI
789 +      }
790 + #endif //is_mpi      
791 +
792 +      molIndex++;
793      }
794    }
795  
# Line 584 | Line 823 | void SimSetup::initFromBass( void ){
823      n_per_extra = (int)ceil( temp1 );
824  
825      if( n_per_extra > 4){
826 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
827 <      exit(8);
826 >      sprintf( painCave.errMsg,
827 >               "SimSetup error. There has been an error in constructing"
828 >               " the non-complete lattice.\n" );
829 >      painCave.isFatal = 1;
830 >      simError();
831      }
832    }
833    else{
# Line 690 | Line 932 | void SimSetup::makeElement( double x, double y, double
932  
933      current_atom = comp_stamps[current_comp]->getAtom( k );
934      if( !current_atom->havePosition() ){
935 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
936 <                << ", atom " << current_atom->getType()
937 <                << " does not have a position specified.\n"
938 <                << "The initialization routine is unable to give a start"
939 <                << " position.\n";
940 <      exit(8);
935 >      sprintf( painCave.errMsg,
936 >               "SimSetup:initFromBass error.\n"
937 >               "\tComponent %s, atom %s does not have a position specified.\n"
938 >               "\tThe initialization routine is unable to give a start"
939 >               " position.\n",
940 >               comp_stamps[current_comp]->getID(),
941 >               current_atom->getType() );
942 >      painCave.isFatal = 1;
943 >      simError();
944      }
945  
946      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );

Diff Legend

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