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 128 by chuckv, Thu Oct 3 21:52:46 2002 UTC vs.
Revision 249 by chuckv, Mon Jan 27 21:28:19 2003 UTC

# Line 6 | Line 6
6   #include "parse_me.h"
7   #include "LRI.hpp"
8   #include "Integrator.hpp"
9 < #include "mpiInterface.h"
9 > #include "simError.h"
10  
11 + #ifdef IS_MPI
12 + #include "mpiBASS.h"
13 + #include "mpiSimulation.hpp"
14 + #include "bassDiag.hpp"
15 + #endif
16 +
17   SimSetup::SimSetup(){
18    stamps = new MakeStamps();
19    globals = new Globals();
20 +  
21 + #ifdef IS_MPI
22 +  strcpy( checkPointMsg, "SimSetup creation successful" );
23 +  MPIcheckPoint();
24 + #endif // IS_MPI
25   }
26  
27   SimSetup::~SimSetup(){
# Line 20 | Line 31 | void SimSetup::parseFile( char* fileName ){
31  
32   void SimSetup::parseFile( char* fileName ){
33  
34 <  inFileName = fileName;
35 <  set_interface_stamps( stamps, globals );
36 < #ifdef MPI
37 <  mpiEventInit();
34 > #ifdef IS_MPI
35 >  if( worldRank == 0 ){
36 > #endif // is_mpi
37 >    
38 >    inFileName = fileName;
39 >    set_interface_stamps( stamps, globals );
40 >    
41 > #ifdef IS_MPI
42 >    mpiEventInit();
43   #endif
44 <  yacc_BASS( fileName );
45 < #ifdef MPI
46 <  throwMPIEvent(NULL);
44 >
45 >    yacc_BASS( fileName );
46 >
47 > #ifdef IS_MPI
48 >    throwMPIEvent(NULL);
49 >  }
50 >  else receiveParse();
51   #endif
52  
53   }
54  
55 < #ifdef MPI
55 > #ifdef IS_MPI
56   void SimSetup::receiveParse(void){
57  
58      set_interface_stamps( stamps, globals );
59      mpiEventInit();
60 +    MPIcheckPoint();
61      mpiEventLoop();
62  
63   }
43 – #endif
64  
65 +
66   void SimSetup::testMe(void){
67    bassDiag* dumpMe = new bassDiag(globals,stamps);
68    dumpMe->dumpStamps();
69    delete dumpMe;
70   }
71 + #endif
72  
73   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 71 | 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 <    std::cerr<< "SimSetup Error. Unrecognized force field -> "
99 <             << force_field << "\n";
100 <    exit(8);
98 >    sprintf( painCave.errMsg,
99 >             "SimSetup Error. Unrecognized force field -> %s\n",
100 >             force_field );
101 >    painCave.isFatal = 1;
102 >    simError();
103    }
104  
105 + #ifdef IS_MPI
106 +  strcpy( checkPointMsg, "ForceField creation successful" );
107 +  MPIcheckPoint();
108 + #endif // is_mpi
109 +
110    // get the components and calculate the tot_nMol and indvidual n_mol
111    the_components = the_globals->getComponents();
112    components_nmol = new int[n_components];
# Line 91 | 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 101 | 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 134 | 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 148 | 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 +  // divide the molecules among processors here.
238 +  
239 +  mpiSim = new mpiSimulation( simnfo );
240 +  
241 +  mpiSim->divideLabor();
242 +
243 +  // set up the local variables
244 +  
245 +  int localMol, allMol;
246 +  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
247 +  
248 +  allMol = 0;
249 +  localMol = 0;
250 +  local_atoms = 0;
251 +  local_bonds = 0;
252 +  local_bends = 0;
253 +  local_torsions = 0;
254 +  for( i=0; i<n_components; i++ ){
255 +
256 +    for( j=0; j<components_nmol[i]; j++ ){
257 +      
258 +      if( mpiSim->getMyMolStart() <= allMol &&
259 +          allMol <= mpiSim->getMyMolEnd() ){
260 +        
261 +        local_atoms +=    comp_stamps[i]->getNAtoms();
262 +        local_bonds +=    comp_stamps[i]->getNBonds();
263 +        local_bends +=    comp_stamps[i]->getNBends();
264 +        local_torsions += comp_stamps[i]->getNTorsions();
265 +        localMol++;
266 +      }      
267 +      allMol++;
268 +    }
269 +  }
270 +  local_SRI = local_bonds + local_bends + local_torsions;
271 +  
272 +
273 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
274 +  
275 +  if( local_atoms != simnfo->n_atoms ){
276 +    sprintf( painCave.errMsg,
277 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
278 +             " localAtom (%d) are note equal.\n",
279 +             simnfo->n_atoms,
280 +             local_atoms );
281 +    painCave.isFatal = 1;
282 +    simError();
283 +  }
284 +
285 +  simnfo->n_bonds = local_bonds;
286 +  simnfo->n_bends = local_bends;
287 +  simnfo->n_torsions = local_torsions;
288 +  simnfo->n_SRI = local_SRI;
289 +  simnfo->n_mol = localMol;
290 +
291 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
292 +  MPIcheckPoint();
293 +  
294 +  
295 + #endif // is_mpi
296 +  
297 +
298    // create the atom and short range interaction arrays
299  
300 <  the_atoms = new Atom*[tot_atoms];
301 <  the_molecules = new Molecule[tot_nmol];
300 >  Atom::createArrays(simnfo->n_atoms);
301 >  the_atoms = new Atom*[simnfo->n_atoms];
302 >  the_molecules = new Molecule[simnfo->n_mol];
303  
304  
305 <  if( tot_SRI ){
306 <    the_sris = new SRI*[tot_SRI];
307 <    the_excludes = new ex_pair[tot_SRI];
305 >  if( simnfo->n_SRI ){
306 >    the_sris = new SRI*[simnfo->n_SRI];
307 >    the_excludes = new ex_pair[simnfo->n_SRI];
308    }
309  
310    // set the arrays into the SimInfo object
# Line 168 | Line 315 | void SimSetup::createSim( void ){
315    simnfo->excludes = the_excludes;
316  
317  
171 –  // initialize the arrays
172 –
173 –  the_ff->setSimInfo( simnfo );
174 –
175 –  makeAtoms();
176 –
177 –  if( tot_bonds ){
178 –    makeBonds();
179 –  }
180 –
181 –  if( tot_bends ){
182 –    makeBends();
183 –  }
184 –
185 –  if( tot_torsions ){
186 –    makeTorsions();
187 –  }
188 –
189 –  //  makeMolecules();
190 –
318    // get some of the tricky things that may still be in the globals
319  
320    if( simnfo->n_dipoles ){
321  
322      if( !the_globals->haveRRF() ){
323 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
324 <      exit(8);
323 >      sprintf( painCave.errMsg,
324 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
325 >      painCave.isFatal = 1;
326 >      simError();
327      }
328      if( !the_globals->haveDielectric() ){
329 <      std::cerr << "SimSetup Error, system has dipoles, but no"
330 <                << " dielectric was set.\n";
331 <      exit(8);
329 >      sprintf( painCave.errMsg,
330 >               "SimSetup Error, system has dipoles, but no"
331 >               " dielectric was set.\n" );
332 >      painCave.isFatal = 1;
333 >      simError();
334      }
335  
336      simnfo->rRF        = the_globals->getRRF();
337      simnfo->dielectric = the_globals->getDielectric();
338    }
339  
340 + #ifdef IS_MPI
341 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
342 +  MPIcheckPoint();
343 + #endif // is_mpi
344 +  
345    if( the_globals->haveBox() ){
346      simnfo->box_x = the_globals->getBox();
347      simnfo->box_y = the_globals->getBox();
# Line 221 | Line 357 | void SimSetup::createSim( void ){
357    }
358    else{
359      if( !the_globals->haveBoxX() ){
360 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
361 <      exit(8);
360 >      sprintf( painCave.errMsg,
361 >               "SimSetup error, no periodic BoxX size given.\n" );
362 >      painCave.isFatal = 1;
363 >      simError();
364      }
365      simnfo->box_x = the_globals->getBoxX();
366  
367      if( !the_globals->haveBoxY() ){
368 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
369 <      exit(8);
368 >      sprintf( painCave.errMsg,
369 >               "SimSetup error, no periodic BoxY size given.\n" );
370 >      painCave.isFatal = 1;
371 >      simError();
372      }
373      simnfo->box_y = the_globals->getBoxY();
374  
375      if( !the_globals->haveBoxZ() ){
376 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
377 <      exit(8);
376 >      sprintf( painCave.errMsg,
377 >               "SimSetup error, no periodic BoxZ size given.\n" );
378 >      painCave.isFatal = 1;
379 >      simError();
380      }
381      simnfo->box_z = the_globals->getBoxZ();
382    }
383  
384 + #ifdef IS_MPI
385 +  strcpy( checkPointMsg, "Box size set up" );
386 +  MPIcheckPoint();
387 + #endif // is_mpi
388  
243 – //   if( the_globals->haveInitialConfig() ){
244 – //        InitializeFromFile* fileInit;
245 – //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
389  
390 < //     fileInit->read_xyz( simnfo ); // default velocities on
390 >  // initialize the arrays
391  
392 < //     delete fileInit;
250 < //   }
251 < //   else{
392 >  the_ff->setSimInfo( simnfo );
393  
394 <    initFromBass();
394 >  makeAtoms();
395  
396 +  if( tot_bonds ){
397 +    makeBonds();
398 +  }
399  
400 < //   }
400 >  if( tot_bends ){
401 >    makeBends();
402 >  }
403  
404 < //   if( the_globals->haveFinalConfig() ){
405 < //     strcpy( simnfo->finalName, the_globals->getFinalConfig() );
406 < //   }
261 < //   else{
262 < //     strcpy( simnfo->finalName, inFileName );
263 < //     char* endTest;
264 < //     int nameLength = strlen( simnfo->finalName );
265 < //     endTest = &(simnfo->finalName[nameLength - 5]);
266 < //     if( !strcmp( endTest, ".bass" ) ){
267 < //       strcpy( endTest, ".eor" );
268 < //     }
269 < //     else if( !strcmp( endTest, ".BASS" ) ){
270 < //       strcpy( endTest, ".eor" );
271 < //     }
272 < //     else{
273 < //       endTest = &(simnfo->finalName[nameLength - 4]);
274 < //       if( !strcmp( endTest, ".bss" ) ){
275 < //      strcpy( endTest, ".eor" );
276 < //       }
277 < //       else if( !strcmp( endTest, ".mdl" ) ){
278 < //      strcpy( endTest, ".eor" );
279 < //       }
280 < //       else{
281 < //      strcat( simnfo->finalName, ".eor" );
282 < //       }
283 < //     }
284 < //   }
404 >  if( tot_torsions ){
405 >    makeTorsions();
406 >  }
407  
286 – //   // make the sample and status out names
408  
288 – //   strcpy( simnfo->sampleName, inFileName );
289 – //   char* endTest;
290 – //   int nameLength = strlen( simnfo->sampleName );
291 – //   endTest = &(simnfo->sampleName[nameLength - 5]);
292 – //   if( !strcmp( endTest, ".bass" ) ){
293 – //     strcpy( endTest, ".dump" );
294 – //   }
295 – //   else if( !strcmp( endTest, ".BASS" ) ){
296 – //     strcpy( endTest, ".dump" );
297 – //   }
298 – //   else{
299 – //     endTest = &(simnfo->sampleName[nameLength - 4]);
300 – //     if( !strcmp( endTest, ".bss" ) ){
301 – //       strcpy( endTest, ".dump" );
302 – //     }
303 – //     else if( !strcmp( endTest, ".mdl" ) ){
304 – //       strcpy( endTest, ".dump" );
305 – //     }
306 – //     else{
307 – //       strcat( simnfo->sampleName, ".dump" );
308 – //     }
309 – //   }
409  
410 < //   strcpy( simnfo->statusName, inFileName );
312 < //   nameLength = strlen( simnfo->statusName );
313 < //   endTest = &(simnfo->statusName[nameLength - 5]);
314 < //   if( !strcmp( endTest, ".bass" ) ){
315 < //     strcpy( endTest, ".stat" );
316 < //   }
317 < //   else if( !strcmp( endTest, ".BASS" ) ){
318 < //     strcpy( endTest, ".stat" );
319 < //   }
320 < //   else{
321 < //     endTest = &(simnfo->statusName[nameLength - 4]);
322 < //     if( !strcmp( endTest, ".bss" ) ){
323 < //       strcpy( endTest, ".stat" );
324 < //     }
325 < //     else if( !strcmp( endTest, ".mdl" ) ){
326 < //       strcpy( endTest, ".stat" );
327 < //     }
328 < //     else{
329 < //       strcat( simnfo->statusName, ".stat" );
330 < //     }
331 < //   }
410 >
411  
412  
413 <  // set the status, sample, and themal kick times
413 > if( the_globals->haveInitialConfig() ){
414 >
415 >     InitializeFromFile* fileInit;
416 > #ifdef IS_MPI // is_mpi
417 >     if( worldRank == 0 ){
418 > #endif //is_mpi
419 >   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
420 > #ifdef IS_MPI
421 >     }else fileInit = new InitializeFromFile( NULL );
422 > #endif
423 >   fileInit->read_xyz( simnfo ); // default velocities on
424  
425 +   delete fileInit;
426 + }
427 + else{
428 +
429 + #ifdef IS_MPI
430 +
431 +  // no init from bass
432 +  
433 +  sprintf( painCave.errMsg,
434 +           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
435 +  painCave.isFatal;
436 +  simError();
437 +  
438 + #else
439 +
440 +  initFromBass();
441 +
442 +
443 + #endif
444 + }
445 +
446 + #ifdef IS_MPI
447 +  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
448 +  MPIcheckPoint();
449 + #endif // is_mpi
450 +
451 +
452 +  
453 +
454 +  
455 +
456 +  
457 + #ifdef IS_MPI
458 +  if( worldRank == 0 ){
459 + #endif // is_mpi
460 +    
461 +    if( the_globals->haveFinalConfig() ){
462 +      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
463 +    }
464 +    else{
465 +      strcpy( simnfo->finalName, inFileName );
466 +      char* endTest;
467 +      int nameLength = strlen( simnfo->finalName );
468 +      endTest = &(simnfo->finalName[nameLength - 5]);
469 +      if( !strcmp( endTest, ".bass" ) ){
470 +        strcpy( endTest, ".eor" );
471 +      }
472 +      else if( !strcmp( endTest, ".BASS" ) ){
473 +        strcpy( endTest, ".eor" );
474 +      }
475 +      else{
476 +        endTest = &(simnfo->finalName[nameLength - 4]);
477 +        if( !strcmp( endTest, ".bss" ) ){
478 +          strcpy( endTest, ".eor" );
479 +        }
480 +        else if( !strcmp( endTest, ".mdl" ) ){
481 +          strcpy( endTest, ".eor" );
482 +        }
483 +        else{
484 +          strcat( simnfo->finalName, ".eor" );
485 +        }
486 +      }
487 +    }
488 +    
489 +    // make the sample and status out names
490 +    
491 +    strcpy( simnfo->sampleName, inFileName );
492 +    char* endTest;
493 +    int nameLength = strlen( simnfo->sampleName );
494 +    endTest = &(simnfo->sampleName[nameLength - 5]);
495 +    if( !strcmp( endTest, ".bass" ) ){
496 +      strcpy( endTest, ".dump" );
497 +    }
498 +    else if( !strcmp( endTest, ".BASS" ) ){
499 +      strcpy( endTest, ".dump" );
500 +    }
501 +    else{
502 +      endTest = &(simnfo->sampleName[nameLength - 4]);
503 +      if( !strcmp( endTest, ".bss" ) ){
504 +        strcpy( endTest, ".dump" );
505 +      }
506 +      else if( !strcmp( endTest, ".mdl" ) ){
507 +        strcpy( endTest, ".dump" );
508 +      }
509 +      else{
510 +        strcat( simnfo->sampleName, ".dump" );
511 +      }
512 +    }
513 +    
514 +    strcpy( simnfo->statusName, inFileName );
515 +    nameLength = strlen( simnfo->statusName );
516 +    endTest = &(simnfo->statusName[nameLength - 5]);
517 +    if( !strcmp( endTest, ".bass" ) ){
518 +      strcpy( endTest, ".stat" );
519 +    }
520 +    else if( !strcmp( endTest, ".BASS" ) ){
521 +      strcpy( endTest, ".stat" );
522 +    }
523 +    else{
524 +      endTest = &(simnfo->statusName[nameLength - 4]);
525 +      if( !strcmp( endTest, ".bss" ) ){
526 +        strcpy( endTest, ".stat" );
527 +      }
528 +      else if( !strcmp( endTest, ".mdl" ) ){
529 +        strcpy( endTest, ".stat" );
530 +      }
531 +      else{
532 +        strcat( simnfo->statusName, ".stat" );
533 +      }
534 +    }
535 +    
536 + #ifdef IS_MPI
537 +  }
538 + #endif // is_mpi
539 +  
540 +  // set the status, sample, and themal kick times
541 +  
542    if( the_globals->haveSampleTime() ){
543      simnfo->sampleTime = the_globals->getSampleTime();
544      simnfo->statusTime = simnfo->sampleTime;
# Line 357 | Line 563 | void SimSetup::createSim( void ){
563    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
564  
565  
566 <  // make the longe range forces and the integrator
566 > //   // make the longe range forces and the integrator
567  
568 <  new AllLong( simnfo );
568 > //   new AllLong( simnfo );
569  
570 <  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
570 >  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
571    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
572    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
573 +  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
574 +
575   }
576  
577   void SimSetup::makeAtoms( void ){
# Line 372 | Line 580 | void SimSetup::makeAtoms( void ){
580    double ux, uy, uz, uSqr, u;
581    AtomStamp* current_atom;
582    DirectionalAtom* dAtom;
583 <  int molIndex, molStart, molEnd, nMemb;
583 >  int molIndex, molStart, molEnd, nMemb, lMolIndex;
584  
585 <
585 >  lMolIndex = 0;
586    molIndex = 0;
587    index = 0;
588    for( i=0; i<n_components; i++ ){
589  
590      for( j=0; j<components_nmol[i]; j++ ){
591  
592 <      molStart = index;
593 <      nMemb = comp_stamps[i]->getNAtoms();
594 <      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
592 > #ifdef IS_MPI
593 >      if( mpiSim->getMyMolStart() <= molIndex &&
594 >          molIndex <= mpiSim->getMyMolEnd() ){
595 > #endif // is_mpi        
596  
597 <        current_atom = comp_stamps[i]->getAtom( k );
598 <        if( current_atom->haveOrientation() ){
599 <
600 <          dAtom = new DirectionalAtom;
601 <          simnfo->n_oriented++;
602 <          the_atoms[index] = dAtom;
603 <
604 <          ux = current_atom->getOrntX();
605 <          uy = current_atom->getOrntY();
606 <          uz = current_atom->getOrntZ();
607 <
608 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
609 <
610 <          u = sqrt( uSqr );
611 <          ux = ux / u;
612 <          uy = uy / u;
613 <          uz = uz / u;
614 <
615 <          dAtom->setSUx( ux );
616 <          dAtom->setSUy( uy );
617 <          dAtom->setSUz( uz );
618 <        }
619 <        else{
620 <          the_atoms[index] = new GeneralAtom;
597 >        molStart = index;
598 >        nMemb = comp_stamps[i]->getNAtoms();
599 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
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          }
632 <        the_atoms[index]->setType( current_atom->getType() );
633 <        the_atoms[index]->setIndex( index );
632 >        
633 >        molEnd = index -1;
634 >        the_molecules[lMolIndex].setNMembers( nMemb );
635 >        the_molecules[lMolIndex].setStartAtom( molStart );
636 >        the_molecules[lMolIndex].setEndAtom( molEnd );
637 >        the_molecules[lMolIndex].setStampID( i );
638 >        lMolIndex++;
639  
640 <        // increment the index and repeat;
417 <        index++;
640 > #ifdef IS_MPI
641        }
642 <
643 <      molEnd = index -1;
421 <      the_molecules[molIndex].setNMembers( nMemb );
422 <      the_molecules[molIndex].setStartAtom( molStart );
423 <      the_molecules[molIndex].setEndAtom( molEnd );
642 > #endif //is_mpi
643 >      
644        molIndex++;
425 –
645      }
646    }
647  
# Line 431 | Line 650 | void SimSetup::makeBonds( void ){
650  
651   void SimSetup::makeBonds( void ){
652  
653 <  int i, j, k, index, offset;
653 >  int i, j, k, index, offset, molIndex;
654    bond_pair* the_bonds;
655    BondStamp* current_bond;
656  
657    the_bonds = new bond_pair[tot_bonds];
658    index = 0;
659    offset = 0;
660 +  molIndex = 0;
661 +
662    for( i=0; i<n_components; i++ ){
663  
664      for( j=0; j<components_nmol[i]; j++ ){
665  
666 <      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
667 <
668 <        current_bond = comp_stamps[i]->getBond( k );
669 <        the_bonds[index].a = current_bond->getA() + offset;
670 <        the_bonds[index].b = current_bond->getB() + offset;
671 <
672 <        the_excludes[index].i = the_bonds[index].a;
673 <        the_excludes[index].j = the_bonds[index].b;
674 <
675 <        // increment the index and repeat;
676 <        index++;
666 > #ifdef IS_MPI
667 >      if( mpiSim->getMyMolStart() <= molIndex &&
668 >          molIndex <= mpiSim->getMyMolEnd() ){
669 > #endif // is_mpi        
670 >        
671 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
672 >          
673 >          current_bond = comp_stamps[i]->getBond( k );
674 >          the_bonds[index].a = current_bond->getA() + offset;
675 >          the_bonds[index].b = current_bond->getB() + offset;
676 >          
677 >          the_excludes[index].i = the_bonds[index].a;
678 >          the_excludes[index].j = the_bonds[index].b;
679 >          
680 >          // increment the index and repeat;
681 >          index++;
682 >        }
683 >        offset += comp_stamps[i]->getNAtoms();
684 >        
685 > #ifdef IS_MPI
686        }
687 <      offset += comp_stamps[i]->getNAtoms();
688 <    }
687 > #endif is_mpi
688 >      
689 >      molIndex++;
690 >    }      
691    }
692  
693    the_ff->initializeBonds( the_bonds );
# Line 463 | Line 695 | void SimSetup::makeBends( void ){
695  
696   void SimSetup::makeBends( void ){
697  
698 <  int i, j, k, index, offset;
698 >  int i, j, k, index, offset, molIndex;
699    bend_set* the_bends;
700    BendStamp* current_bend;
701  
702    the_bends = new bend_set[tot_bends];
703    index = 0;
704    offset = 0;
705 +  molIndex = 0;
706    for( i=0; i<n_components; i++ ){
707  
708      for( j=0; j<components_nmol[i]; j++ ){
709  
710 <      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
710 > #ifdef IS_MPI
711 >      if( mpiSim->getMyMolStart() <= molIndex &&
712 >          molIndex <= mpiSim->getMyMolEnd() ){
713 > #endif // is_mpi        
714  
715 <        current_bend = comp_stamps[i]->getBend( k );
716 <        the_bends[index].a = current_bend->getA() + offset;
717 <        the_bends[index].b = current_bend->getB() + offset;
718 <        the_bends[index].c = current_bend->getC() + offset;
719 <
720 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
721 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
722 <
723 <        // increment the index and repeat;
724 <        index++;
715 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
716 >          
717 >          current_bend = comp_stamps[i]->getBend( k );
718 >          the_bends[index].a = current_bend->getA() + offset;
719 >          the_bends[index].b = current_bend->getB() + offset;
720 >          the_bends[index].c = current_bend->getC() + offset;
721 >          
722 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
723 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
724 >          
725 >          // increment the index and repeat;
726 >          index++;
727 >        }
728 >        offset += comp_stamps[i]->getNAtoms();
729 >        
730 > #ifdef IS_MPI
731        }
732 <      offset += comp_stamps[i]->getNAtoms();
732 > #endif //is_mpi
733 >
734 >      molIndex++;
735      }
736    }
737  
# Line 496 | Line 740 | void SimSetup::makeTorsions( void ){
740  
741   void SimSetup::makeTorsions( void ){
742  
743 <  int i, j, k, index, offset;
743 >  int i, j, k, index, offset, molIndex;
744    torsion_set* the_torsions;
745    TorsionStamp* current_torsion;
746  
747    the_torsions = new torsion_set[tot_torsions];
748    index = 0;
749    offset = 0;
750 +  molIndex = 0;
751    for( i=0; i<n_components; i++ ){
752  
753      for( j=0; j<components_nmol[i]; j++ ){
754  
755 + #ifdef IS_MPI
756 +      if( mpiSim->getMyMolStart() <= molIndex &&
757 +          molIndex <= mpiSim->getMyMolEnd() ){
758 + #endif // is_mpi        
759 +
760        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
761  
762          current_torsion = comp_stamps[i]->getTorsion( k );
# Line 522 | Line 772 | void SimSetup::makeTorsions( void ){
772          index++;
773        }
774        offset += comp_stamps[i]->getNAtoms();
775 +
776 + #ifdef IS_MPI
777 +      }
778 + #endif //is_mpi      
779 +
780 +      molIndex++;
781      }
782    }
783  
# Line 555 | Line 811 | void SimSetup::initFromBass( void ){
811      n_per_extra = (int)ceil( temp1 );
812  
813      if( n_per_extra > 4){
814 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
815 <      exit(8);
814 >      sprintf( painCave.errMsg,
815 >               "SimSetup error. There has been an error in constructing"
816 >               " the non-complete lattice.\n" );
817 >      painCave.isFatal = 1;
818 >      simError();
819      }
820    }
821    else{
# Line 661 | Line 920 | void SimSetup::makeElement( double x, double y, double
920  
921      current_atom = comp_stamps[current_comp]->getAtom( k );
922      if( !current_atom->havePosition() ){
923 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
924 <                << ", atom " << current_atom->getType()
925 <                << " does not have a position specified.\n"
926 <                << "The initialization routine is unable to give a start"
927 <                << " position.\n";
928 <      exit(8);
923 >      sprintf( painCave.errMsg,
924 >               "SimSetup:initFromBass error.\n"
925 >               "\tComponent %s, atom %s does not have a position specified.\n"
926 >               "\tThe initialization routine is unable to give a start"
927 >               " position.\n",
928 >               comp_stamps[current_comp]->getID(),
929 >               current_atom->getType() );
930 >      painCave.isFatal = 1;
931 >      simError();
932      }
933  
934      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );

Diff Legend

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