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

Comparing:
branches/mmeineke/mdtools/interface_implementation/SimSetup.cpp (file contents), Revision 10 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
trunk/mdtools/interface_implementation/SimSetup.cpp (file contents), Revision 194 by chuckv, Wed Dec 4 21:19:38 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;
# 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  
165    for( i=0; i<n_components; i++ ){
166  
167 <    comp_stamps[i] =
167 >    comp_stamps[i] =
168        the_stamps->getMolecule( the_components[i]->getType() );
169    }
170  
104 –  
171  
172 +
173    // caclulate the number of atoms, bonds, bends and torsions
174  
175    tot_atoms = 0;
# Line 110 | Line 177 | void SimSetup::createSim( void ){
177    tot_bends = 0;
178    tot_torsions = 0;
179    for( i=0; i<n_components; i++ ){
180 <    
180 >
181      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
182      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
183      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
184      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
185    }
186 <  
186 >
187    tot_SRI = tot_bonds + tot_bends + tot_torsions;
188 <  
188 >
189    simnfo->n_atoms = tot_atoms;
190    simnfo->n_bonds = tot_bonds;
191    simnfo->n_bends = tot_bends;
192    simnfo->n_torsions = tot_torsions;
193    simnfo->n_SRI = tot_SRI;
194  
195 +  // divide the molecules among processors here.
196 +
197 +
198    // create the atom and short range interaction arrays
199 <  
199 >
200 >  Atom::createArrays(tot_atoms);
201    the_atoms = new Atom*[tot_atoms];
202 <  //  the_molecules = new Molecule[tot_nmol];
203 <  
204 <  
202 >  the_molecules = new Molecule[tot_nmol];
203 >
204 >
205    if( tot_SRI ){
206      the_sris = new SRI*[tot_SRI];
207      the_excludes = new ex_pair[tot_SRI];
# Line 143 | Line 214 | void SimSetup::createSim( void ){
214    simnfo->n_exclude = tot_SRI;
215    simnfo->excludes = the_excludes;
216  
217 +
218    // initialize the arrays
219 <  
219 >
220    the_ff->setSimInfo( simnfo );
221 <    
221 >
222    makeAtoms();
223  
224    if( tot_bonds ){
# Line 168 | Line 240 | void SimSetup::createSim( void ){
240    if( simnfo->n_dipoles ){
241  
242      if( !the_globals->haveRRF() ){
243 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
244 <      exit(8);
243 >      sprintf( painCave.errMsg,
244 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
245 >      painCave.isFatal = 1;
246 >      simError();
247      }
248      if( !the_globals->haveDielectric() ){
249 <      std::cerr << "SimSetup Error, system has dipoles, but no"
250 <                << " dielectric was set.\n";
251 <      exit(8);
249 >      sprintf( painCave.errMsg,
250 >               "SimSetup Error, system has dipoles, but no"
251 >               " dielectric was set.\n" );
252 >      painCave.isFatal = 1;
253 >      simError();
254      }
255  
256      simnfo->rRF        = the_globals->getRRF();
257      simnfo->dielectric = the_globals->getDielectric();
258    }
259  
260 + #ifdef IS_MPI
261 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
262 +  MPIcheckPoint();
263 + #endif // is_mpi
264 +  
265    if( the_globals->haveBox() ){
266      simnfo->box_x = the_globals->getBox();
267      simnfo->box_y = the_globals->getBox();
268      simnfo->box_z = the_globals->getBox();
269    }
270    else if( the_globals->haveDensity() ){
271 <    
271 >
272      double vol;
273      vol = (double)tot_nmol / the_globals->getDensity();
274      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 196 | Line 277 | void SimSetup::createSim( void ){
277    }
278    else{
279      if( !the_globals->haveBoxX() ){
280 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
281 <      exit(8);
280 >      sprintf( painCave.errMsg,
281 >               "SimSetup error, no periodic BoxX size given.\n" );
282 >      painCave.isFatal = 1;
283 >      simError();
284      }
285      simnfo->box_x = the_globals->getBoxX();
286  
287      if( !the_globals->haveBoxY() ){
288 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
289 <      exit(8);
288 >      sprintf( painCave.errMsg,
289 >               "SimSetup error, no periodic BoxY size given.\n" );
290 >      painCave.isFatal = 1;
291 >      simError();
292      }
293      simnfo->box_y = the_globals->getBoxY();
294  
295      if( !the_globals->haveBoxZ() ){
296 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
297 <      exit(8);
296 >      sprintf( painCave.errMsg,
297 >               "SimSetup error, no periodic BoxZ size given.\n" );
298 >      painCave.isFatal = 1;
299 >      simError();
300      }
301      simnfo->box_z = the_globals->getBoxZ();
302    }
216 –    
217 –  if( the_globals->haveInitialConfig() ){
218 –    InitializeFromFile* fileInit;
219 –    fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
220 –    
221 –    fileInit->read_xyz( simnfo ); // default velocities on
303  
304 <    delete fileInit;    
305 <  }
306 <  else{
307 <    initFromBass();
227 <  }
304 > #ifdef IS_MPI
305 >  strcpy( checkPointMsg, "Box size set up" );
306 >  MPIcheckPoint();
307 > #endif // is_mpi
308  
309 <  if( the_globals->haveFinalConfig() ){
310 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
311 <  }
312 <  else{
313 <    strcpy( simnfo->finalName, inFileName );
314 <    char* endTest;
315 <    int nameLength = strlen( simnfo->finalName );
316 <    endTest = &(simnfo->finalName[nameLength - 5]);
317 <    if( !strcmp( endTest, ".bass" ) ){
318 <      strcpy( endTest, ".eor" );
309 >
310 >
311 > //   if( the_globals->haveInitialConfig() ){
312 > //        InitializeFromFile* fileInit;
313 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
314 >
315 > //     fileInit->read_xyz( simnfo ); // default velocities on
316 >
317 > //     delete fileInit;
318 > //   }
319 > //   else{
320 >
321 >  initFromBass();
322 >
323 > #ifdef IS_MPI
324 >  strcpy( checkPointMsg, "initFromBass successfully created the lattice" );
325 >  MPIcheckPoint();
326 > #endif // is_mpi
327 >
328 >
329 >  
330 >
331 >  
332 >  //   }
333 >  
334 > #ifdef IS_MPI
335 >  if( worldRank == 0 ){
336 > #endif // is_mpi
337 >    
338 >    if( the_globals->haveFinalConfig() ){
339 >      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
340      }
240 –    else if( !strcmp( endTest, ".BASS" ) ){
241 –      strcpy( endTest, ".eor" );
242 –    }
341      else{
342 <      endTest = &(simnfo->finalName[nameLength - 4]);
343 <      if( !strcmp( endTest, ".bss" ) ){
342 >      strcpy( simnfo->finalName, inFileName );
343 >      char* endTest;
344 >      int nameLength = strlen( simnfo->finalName );
345 >      endTest = &(simnfo->finalName[nameLength - 5]);
346 >      if( !strcmp( endTest, ".bass" ) ){
347          strcpy( endTest, ".eor" );
348        }
349 <      else if( !strcmp( endTest, ".mdl" ) ){
349 >      else if( !strcmp( endTest, ".BASS" ) ){
350          strcpy( endTest, ".eor" );
351        }
352        else{
353 <        strcat( simnfo->finalName, ".eor" );
353 >        endTest = &(simnfo->finalName[nameLength - 4]);
354 >        if( !strcmp( endTest, ".bss" ) ){
355 >          strcpy( endTest, ".eor" );
356 >        }
357 >        else if( !strcmp( endTest, ".mdl" ) ){
358 >          strcpy( endTest, ".eor" );
359 >        }
360 >        else{
361 >          strcat( simnfo->finalName, ".eor" );
362 >        }
363        }
364      }
365 <  }
365 >    
366 >    // make the sample and status out names
367      
368 <  // make the sample and status out names
369 <
370 <  strcpy( simnfo->sampleName, inFileName );
371 <  char* endTest;
372 <  int nameLength = strlen( simnfo->sampleName );
262 <  endTest = &(simnfo->sampleName[nameLength - 5]);
263 <  if( !strcmp( endTest, ".bass" ) ){
264 <    strcpy( endTest, ".dump" );
265 <  }
266 <  else if( !strcmp( endTest, ".BASS" ) ){
267 <    strcpy( endTest, ".dump" );
268 <  }
269 <  else{
270 <    endTest = &(simnfo->sampleName[nameLength - 4]);
271 <    if( !strcmp( endTest, ".bss" ) ){
368 >    strcpy( simnfo->sampleName, inFileName );
369 >    char* endTest;
370 >    int nameLength = strlen( simnfo->sampleName );
371 >    endTest = &(simnfo->sampleName[nameLength - 5]);
372 >    if( !strcmp( endTest, ".bass" ) ){
373        strcpy( endTest, ".dump" );
374      }
375 <    else if( !strcmp( endTest, ".mdl" ) ){
375 >    else if( !strcmp( endTest, ".BASS" ) ){
376        strcpy( endTest, ".dump" );
377      }
378      else{
379 <      strcat( simnfo->sampleName, ".dump" );
379 >      endTest = &(simnfo->sampleName[nameLength - 4]);
380 >      if( !strcmp( endTest, ".bss" ) ){
381 >        strcpy( endTest, ".dump" );
382 >      }
383 >      else if( !strcmp( endTest, ".mdl" ) ){
384 >        strcpy( endTest, ".dump" );
385 >      }
386 >      else{
387 >        strcat( simnfo->sampleName, ".dump" );
388 >      }
389      }
390 <  }
391 <  
392 <  strcpy( simnfo->statusName, inFileName );
393 <  nameLength = strlen( simnfo->statusName );
394 <  endTest = &(simnfo->statusName[nameLength - 5]);
285 <  if( !strcmp( endTest, ".bass" ) ){
286 <    strcpy( endTest, ".stat" );
287 <  }
288 <  else if( !strcmp( endTest, ".BASS" ) ){
289 <    strcpy( endTest, ".stat" );
290 <  }
291 <  else{
292 <    endTest = &(simnfo->statusName[nameLength - 4]);
293 <    if( !strcmp( endTest, ".bss" ) ){
390 >    
391 >    strcpy( simnfo->statusName, inFileName );
392 >    nameLength = strlen( simnfo->statusName );
393 >    endTest = &(simnfo->statusName[nameLength - 5]);
394 >    if( !strcmp( endTest, ".bass" ) ){
395        strcpy( endTest, ".stat" );
396      }
397 <    else if( !strcmp( endTest, ".mdl" ) ){
397 >    else if( !strcmp( endTest, ".BASS" ) ){
398        strcpy( endTest, ".stat" );
399      }
400      else{
401 <      strcat( simnfo->statusName, ".stat" );
401 >      endTest = &(simnfo->statusName[nameLength - 4]);
402 >      if( !strcmp( endTest, ".bss" ) ){
403 >        strcpy( endTest, ".stat" );
404 >      }
405 >      else if( !strcmp( endTest, ".mdl" ) ){
406 >        strcpy( endTest, ".stat" );
407 >      }
408 >      else{
409 >        strcat( simnfo->statusName, ".stat" );
410 >      }
411      }
412 +    
413 + #ifdef IS_MPI
414    }
415 + #endif // is_mpi
416    
417    // set the status, sample, and themal kick times
418 <
418 >  
419    if( the_globals->haveSampleTime() ){
420 <    simnfo->sampleTime = the_globals->getSampleTime();
420 >    simnfo->sampleTime = the_globals->getSampleTime();
421      simnfo->statusTime = simnfo->sampleTime;
422      simnfo->thermalTime = simnfo->sampleTime;
423    }
424    else{
425 <    simnfo->sampleTime = the_globals->getRunTime();
425 >    simnfo->sampleTime = the_globals->getRunTime();
426      simnfo->statusTime = simnfo->sampleTime;
427      simnfo->thermalTime = simnfo->sampleTime;
428    }
# Line 325 | Line 438 | void SimSetup::createSim( void ){
438    // check for the temperature set flag
439  
440    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
441 <  
442 <    
441 >
442 >
443    // make the longe range forces and the integrator
444 <  
444 >
445    new AllLong( simnfo );
446 <      
446 >
447    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
448    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
449    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
450   }
451  
452   void SimSetup::makeAtoms( void ){
453 <  
453 >
454    int i, j, k, index;
455    double ux, uy, uz, uSqr, u;
456    AtomStamp* current_atom;
457    DirectionalAtom* dAtom;
458 +  int molIndex, molStart, molEnd, nMemb;
459  
460 +
461 +  molIndex = 0;
462    index = 0;
463    for( i=0; i<n_components; i++ ){
464 <    
464 >
465      for( j=0; j<components_nmol[i]; j++ ){
466 <      
466 >
467 >      molStart = index;
468 >      nMemb = comp_stamps[i]->getNAtoms();
469        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
470 <        
470 >
471          current_atom = comp_stamps[i]->getAtom( k );
472 <        if( current_atom->haveOrientation() ){
472 >        if( current_atom->haveOrientation() ){
473  
474 <          dAtom = new DirectionalAtom;
474 >          dAtom = new DirectionalAtom(index);
475            simnfo->n_oriented++;
476            the_atoms[index] = dAtom;
477 <      
477 >
478            ux = current_atom->getOrntX();
479            uy = current_atom->getOrntY();
480            uz = current_atom->getOrntZ();
481 <          
481 >
482            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
483 <          
483 >
484            u = sqrt( uSqr );
485            ux = ux / u;
486            uy = uy / u;
487            uz = uz / u;
488 <          
488 >
489            dAtom->setSUx( ux );
490            dAtom->setSUy( uy );
491            dAtom->setSUz( uz );
492          }
493          else{
494 <          the_atoms[index] = new GeneralAtom;
494 >          the_atoms[index] = new GeneralAtom(index);
495          }
496          the_atoms[index]->setType( current_atom->getType() );
497          the_atoms[index]->setIndex( index );
498 <        
498 >
499          // increment the index and repeat;
500          index++;
501        }
502 +
503 +      molEnd = index -1;
504 +      the_molecules[molIndex].setNMembers( nMemb );
505 +      the_molecules[molIndex].setStartAtom( molStart );
506 +      the_molecules[molIndex].setEndAtom( molEnd );
507 +      molIndex++;
508 +
509      }
510    }
511 <  
511 >
512    the_ff->initializeAtoms();
513   }
514  
# Line 397 | Line 522 | void SimSetup::makeBonds( void ){
522    index = 0;
523    offset = 0;
524    for( i=0; i<n_components; i++ ){
525 <    
525 >
526      for( j=0; j<components_nmol[i]; j++ ){
527 <      
527 >
528        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
529 <        
529 >
530          current_bond = comp_stamps[i]->getBond( k );
531          the_bonds[index].a = current_bond->getA() + offset;
532          the_bonds[index].b = current_bond->getB() + offset;
# Line 415 | Line 540 | void SimSetup::makeBonds( void ){
540        offset += comp_stamps[i]->getNAtoms();
541      }
542    }
543 <  
543 >
544    the_ff->initializeBonds( the_bonds );
545   }
546  
# Line 429 | Line 554 | void SimSetup::makeBends( void ){
554    index = 0;
555    offset = 0;
556    for( i=0; i<n_components; i++ ){
557 <    
557 >
558      for( j=0; j<components_nmol[i]; j++ ){
559 <      
559 >
560        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
561 <        
561 >
562          current_bend = comp_stamps[i]->getBend( k );
563          the_bends[index].a = current_bend->getA() + offset;
564          the_bends[index].b = current_bend->getB() + offset;
# Line 448 | Line 573 | void SimSetup::makeBends( void ){
573        offset += comp_stamps[i]->getNAtoms();
574      }
575    }
576 <  
576 >
577    the_ff->initializeBends( the_bends );
578   }
579  
# Line 462 | Line 587 | void SimSetup::makeTorsions( void ){
587    index = 0;
588    offset = 0;
589    for( i=0; i<n_components; i++ ){
590 <    
590 >
591      for( j=0; j<components_nmol[i]; j++ ){
592 <      
592 >
593        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
594 <        
594 >
595          current_torsion = comp_stamps[i]->getTorsion( k );
596          the_torsions[index].a = current_torsion->getA() + offset;
597          the_torsions[index].b = current_torsion->getB() + offset;
# Line 482 | Line 607 | void SimSetup::makeTorsions( void ){
607        offset += comp_stamps[i]->getNAtoms();
608      }
609    }
610 <  
610 >
611    the_ff->initializeTorsions( the_torsions );
612   }
613  
489 – void SimSetup::makeMolecules( void ){
490 –
491 –  //empy for now
492 – }
493 –
614   void SimSetup::initFromBass( void ){
615  
616    int i, j, k;
# Line 518 | Line 638 | void SimSetup::initFromBass( void ){
638      n_per_extra = (int)ceil( temp1 );
639  
640      if( n_per_extra > 4){
641 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
642 <      exit(8);
641 >      sprintf( painCave.errMsg,
642 >               "SimSetup error. There has been an error in constructing"
643 >               " the non-complete lattice.\n" );
644 >      painCave.isFatal = 1;
645 >      simError();
646      }
647    }
648    else{
# Line 528 | Line 651 | void SimSetup::initFromBass( void ){
651      celly = simnfo->box_y / temp3;
652      cellz = simnfo->box_z / temp3;
653    }
654 <  
654 >
655    current_mol = 0;
656    current_comp_mol = 0;
657    current_comp = 0;
658    current_atom_ndx = 0;
659 <  
659 >
660    for( i=0; i < n_cells ; i++ ){
661      for( j=0; j < n_cells; j++ ){
662        for( k=0; k < n_cells; k++ ){
663 <        
663 >
664          makeElement( i * cellx,
665                       j * celly,
666                       k * cellz );
667 <        
667 >
668          makeElement( i * cellx + 0.5 * cellx,
669                       j * celly + 0.5 * celly,
670                       k * cellz );
671 <        
671 >
672          makeElement( i * cellx,
673                       j * celly + 0.5 * celly,
674                       k * cellz + 0.5 * cellz );
675 <        
675 >
676          makeElement( i * cellx + 0.5 * cellx,
677                       j * celly,
678                       k * cellz + 0.5 * cellz );
# Line 559 | Line 682 | void SimSetup::initFromBass( void ){
682  
683    if( have_extra ){
684      done = 0;
685 <    
685 >
686      int start_ndx;
687      for( i=0; i < (n_cells+1) && !done; i++ ){
688        for( j=0; j < (n_cells+1) && !done; j++ ){
689 <        
689 >
690          if( i < n_cells ){
691 <          
691 >
692            if( j < n_cells ){
693              start_ndx = n_cells;
694            }
695            else start_ndx = 0;
696          }
697          else start_ndx = 0;
698 <        
698 >
699          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
700 <          
700 >
701            makeElement( i * cellx,
702                         j * celly,
703                         k * cellz );
704            done = ( current_mol >= tot_nmol );
705 <          
705 >
706            if( !done && n_per_extra > 1 ){
707              makeElement( i * cellx + 0.5 * cellx,
708                           j * celly + 0.5 * celly,
709                           k * cellz );
710              done = ( current_mol >= tot_nmol );
711            }
712 <          
712 >
713            if( !done && n_per_extra > 2){
714              makeElement( i * cellx,
715                           j * celly + 0.5 * celly,
716                           k * cellz + 0.5 * cellz );
717              done = ( current_mol >= tot_nmol );
718            }
719 <          
719 >
720            if( !done && n_per_extra > 3){
721              makeElement( i * cellx + 0.5 * cellx,
722                           j * celly,
# Line 604 | Line 727 | void SimSetup::initFromBass( void ){
727        }
728      }
729    }
730 <  
731 <  
730 >
731 >
732    for( i=0; i<simnfo->n_atoms; i++ ){
733      simnfo->atoms[i]->set_vx( 0.0 );
734      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 744 | void SimSetup::makeElement( double x, double y, double
744    double rotMat[3][3];
745  
746    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
747 <    
747 >
748      current_atom = comp_stamps[current_comp]->getAtom( k );
749      if( !current_atom->havePosition() ){
750 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
751 <                << ", atom " << current_atom->getType()
752 <                << " does not have a position specified.\n"
753 <                << "The initialization routine is unable to give a start"
754 <                << " position.\n";
755 <      exit(8);
750 >      sprintf( painCave.errMsg,
751 >               "SimSetup:initFromBass error.\n"
752 >               "\tComponent %s, atom %s does not have a position specified.\n"
753 >               "\tThe initialization routine is unable to give a start"
754 >               " position.\n",
755 >               comp_stamps[current_comp]->getID(),
756 >               current_atom->getType() );
757 >      painCave.isFatal = 1;
758 >      simError();
759      }
760 <    
760 >
761      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
762      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
763      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
764 <    
764 >
765      if( the_atoms[current_atom_ndx]->isDirectional() ){
766 <      
766 >
767        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
768 <      
768 >
769        rotMat[0][0] = 1.0;
770        rotMat[0][1] = 0.0;
771        rotMat[0][2] = 0.0;
# Line 657 | Line 783 | void SimSetup::makeElement( double x, double y, double
783  
784      current_atom_ndx++;
785    }
786 <  
786 >
787    current_mol++;
788    current_comp_mol++;
789  
790    if( current_comp_mol >= components_nmol[current_comp] ){
791 <    
791 >
792      current_comp_mol = 0;
793      current_comp++;
794    }

Diff Legend

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