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

Diff Legend

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