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 195 by chuckv, Thu Dec 5 18:53:40 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 +  simnfo->n_mol = tot_nmol;
195  
128 –  // create the atom and short range interaction arrays
196    
197 <  the_atoms = new Atom*[tot_atoms];
198 <  the_molecules = new Molecule[tot_nmol];
197 > #ifdef IS_MPI
198 >
199 >  // divide the molecules among processors here.
200    
201 +  new mpiSimulation( simnfo );
202    
203 <  if( tot_SRI ){
204 <    the_sris = new SRI*[tot_SRI];
205 <    the_excludes = new ex_pair[tot_SRI];
203 >  simnfo->mpiSim->divideLabor( n_components, comp_stamps, components_nmol );
204 >
205 > #endif // is_mpi
206 >  
207 >
208 >  // create the atom and short range interaction arrays
209 >
210 >  Atom::createArrays(simnfo->n_atoms);
211 >  the_atoms = new Atom*[simnfo->n_atoms];
212 >  the_molecules = new Molecule[simnfo->n_mol];
213 >
214 >
215 >  if( simnfo->n_SRI ){
216 >    the_sris = new SRI*[simnfo->n_SRI];
217 >    the_excludes = new ex_pair[simnfo->n_SRI];
218    }
219  
220    // set the arrays into the SimInfo object
# Line 142 | Line 223 | void SimSetup::createSim( void ){
223    simnfo->sr_interactions = the_sris;
224    simnfo->n_exclude = tot_SRI;
225    simnfo->excludes = the_excludes;
145 –  
226  
227 +
228    // initialize the arrays
229 <  
229 >
230    the_ff->setSimInfo( simnfo );
231 <    
231 >
232    makeAtoms();
233  
234    if( tot_bonds ){
# Line 162 | Line 243 | void SimSetup::createSim( void ){
243      makeTorsions();
244    }
245  
165 –  //  makeMolecules();
246  
247    // get some of the tricky things that may still be in the globals
248  
249    if( simnfo->n_dipoles ){
250  
251      if( !the_globals->haveRRF() ){
252 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
253 <      exit(8);
252 >      sprintf( painCave.errMsg,
253 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
254 >      painCave.isFatal = 1;
255 >      simError();
256      }
257      if( !the_globals->haveDielectric() ){
258 <      std::cerr << "SimSetup Error, system has dipoles, but no"
259 <                << " dielectric was set.\n";
260 <      exit(8);
258 >      sprintf( painCave.errMsg,
259 >               "SimSetup Error, system has dipoles, but no"
260 >               " dielectric was set.\n" );
261 >      painCave.isFatal = 1;
262 >      simError();
263      }
264  
265      simnfo->rRF        = the_globals->getRRF();
266      simnfo->dielectric = the_globals->getDielectric();
267    }
268  
269 + #ifdef IS_MPI
270 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
271 +  MPIcheckPoint();
272 + #endif // is_mpi
273 +  
274    if( the_globals->haveBox() ){
275      simnfo->box_x = the_globals->getBox();
276      simnfo->box_y = the_globals->getBox();
277      simnfo->box_z = the_globals->getBox();
278    }
279    else if( the_globals->haveDensity() ){
280 <    
280 >
281      double vol;
282      vol = (double)tot_nmol / the_globals->getDensity();
283      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 197 | Line 286 | void SimSetup::createSim( void ){
286    }
287    else{
288      if( !the_globals->haveBoxX() ){
289 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
290 <      exit(8);
289 >      sprintf( painCave.errMsg,
290 >               "SimSetup error, no periodic BoxX size given.\n" );
291 >      painCave.isFatal = 1;
292 >      simError();
293      }
294      simnfo->box_x = the_globals->getBoxX();
295  
296      if( !the_globals->haveBoxY() ){
297 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
298 <      exit(8);
297 >      sprintf( painCave.errMsg,
298 >               "SimSetup error, no periodic BoxY size given.\n" );
299 >      painCave.isFatal = 1;
300 >      simError();
301      }
302      simnfo->box_y = the_globals->getBoxY();
303  
304      if( !the_globals->haveBoxZ() ){
305 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
306 <      exit(8);
305 >      sprintf( painCave.errMsg,
306 >               "SimSetup error, no periodic BoxZ size given.\n" );
307 >      painCave.isFatal = 1;
308 >      simError();
309      }
310      simnfo->box_z = the_globals->getBoxZ();
311    }
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
312  
313 <    delete fileInit;    
314 <  }
315 <  else{
316 <    initFromBass();
228 <  }
313 > #ifdef IS_MPI
314 >  strcpy( checkPointMsg, "Box size set up" );
315 >  MPIcheckPoint();
316 > #endif // is_mpi
317  
318 <  if( the_globals->haveFinalConfig() ){
319 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
320 <  }
321 <  else{
322 <    strcpy( simnfo->finalName, inFileName );
323 <    char* endTest;
324 <    int nameLength = strlen( simnfo->finalName );
325 <    endTest = &(simnfo->finalName[nameLength - 5]);
326 <    if( !strcmp( endTest, ".bass" ) ){
327 <      strcpy( endTest, ".eor" );
328 <    }
329 <    else if( !strcmp( endTest, ".BASS" ) ){
330 <      strcpy( endTest, ".eor" );
318 >
319 >
320 > //   if( the_globals->haveInitialConfig() ){
321 > //        InitializeFromFile* fileInit;
322 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
323 >
324 > //     fileInit->read_xyz( simnfo ); // default velocities on
325 >
326 > //     delete fileInit;
327 > //   }
328 > //   else{
329 >
330 > #ifdef IS_MPI
331 >
332 >  // no init from bass
333 >  
334 >  sprintf( painCave.errMsg,
335 >           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
336 >  painCave.isFatal;
337 >  simError();
338 >  
339 > #else
340 >
341 >  initFromBass();
342 >
343 > #endif // is_mpi
344 >
345 > #ifdef IS_MPI
346 >  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
347 >  MPIcheckPoint();
348 > #endif // is_mpi
349 >
350 >
351 >  
352 >
353 >  
354 >  //   }
355 >  
356 > #ifdef IS_MPI
357 >  if( worldRank == 0 ){
358 > #endif // is_mpi
359 >    
360 >    if( the_globals->haveFinalConfig() ){
361 >      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
362      }
363      else{
364 <      endTest = &(simnfo->finalName[nameLength - 4]);
365 <      if( !strcmp( endTest, ".bss" ) ){
364 >      strcpy( simnfo->finalName, inFileName );
365 >      char* endTest;
366 >      int nameLength = strlen( simnfo->finalName );
367 >      endTest = &(simnfo->finalName[nameLength - 5]);
368 >      if( !strcmp( endTest, ".bass" ) ){
369          strcpy( endTest, ".eor" );
370        }
371 <      else if( !strcmp( endTest, ".mdl" ) ){
371 >      else if( !strcmp( endTest, ".BASS" ) ){
372          strcpy( endTest, ".eor" );
373        }
374        else{
375 <        strcat( simnfo->finalName, ".eor" );
375 >        endTest = &(simnfo->finalName[nameLength - 4]);
376 >        if( !strcmp( endTest, ".bss" ) ){
377 >          strcpy( endTest, ".eor" );
378 >        }
379 >        else if( !strcmp( endTest, ".mdl" ) ){
380 >          strcpy( endTest, ".eor" );
381 >        }
382 >        else{
383 >          strcat( simnfo->finalName, ".eor" );
384 >        }
385        }
386      }
256 –  }
387      
388 <  // make the sample and status out names
389 <
390 <  strcpy( simnfo->sampleName, inFileName );
391 <  char* endTest;
392 <  int nameLength = strlen( simnfo->sampleName );
393 <  endTest = &(simnfo->sampleName[nameLength - 5]);
394 <  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" ) ){
388 >    // make the sample and status out names
389 >    
390 >    strcpy( simnfo->sampleName, inFileName );
391 >    char* endTest;
392 >    int nameLength = strlen( simnfo->sampleName );
393 >    endTest = &(simnfo->sampleName[nameLength - 5]);
394 >    if( !strcmp( endTest, ".bass" ) ){
395        strcpy( endTest, ".dump" );
396      }
397 <    else if( !strcmp( endTest, ".mdl" ) ){
397 >    else if( !strcmp( endTest, ".BASS" ) ){
398        strcpy( endTest, ".dump" );
399      }
400      else{
401 <      strcat( simnfo->sampleName, ".dump" );
401 >      endTest = &(simnfo->sampleName[nameLength - 4]);
402 >      if( !strcmp( endTest, ".bss" ) ){
403 >        strcpy( endTest, ".dump" );
404 >      }
405 >      else if( !strcmp( endTest, ".mdl" ) ){
406 >        strcpy( endTest, ".dump" );
407 >      }
408 >      else{
409 >        strcat( simnfo->sampleName, ".dump" );
410 >      }
411      }
412 <  }
413 <  
414 <  strcpy( simnfo->statusName, inFileName );
415 <  nameLength = strlen( simnfo->statusName );
416 <  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" ) ){
412 >    
413 >    strcpy( simnfo->statusName, inFileName );
414 >    nameLength = strlen( simnfo->statusName );
415 >    endTest = &(simnfo->statusName[nameLength - 5]);
416 >    if( !strcmp( endTest, ".bass" ) ){
417        strcpy( endTest, ".stat" );
418      }
419 <    else if( !strcmp( endTest, ".mdl" ) ){
419 >    else if( !strcmp( endTest, ".BASS" ) ){
420        strcpy( endTest, ".stat" );
421      }
422      else{
423 <      strcat( simnfo->statusName, ".stat" );
423 >      endTest = &(simnfo->statusName[nameLength - 4]);
424 >      if( !strcmp( endTest, ".bss" ) ){
425 >        strcpy( endTest, ".stat" );
426 >      }
427 >      else if( !strcmp( endTest, ".mdl" ) ){
428 >        strcpy( endTest, ".stat" );
429 >      }
430 >      else{
431 >        strcat( simnfo->statusName, ".stat" );
432 >      }
433      }
434 +    
435 + #ifdef IS_MPI
436    }
437 + #endif // is_mpi
438    
439    // set the status, sample, and themal kick times
440 <
440 >  
441    if( the_globals->haveSampleTime() ){
442 <    simnfo->sampleTime = the_globals->getSampleTime();
442 >    simnfo->sampleTime = the_globals->getSampleTime();
443      simnfo->statusTime = simnfo->sampleTime;
444      simnfo->thermalTime = simnfo->sampleTime;
445    }
446    else{
447 <    simnfo->sampleTime = the_globals->getRunTime();
447 >    simnfo->sampleTime = the_globals->getRunTime();
448      simnfo->statusTime = simnfo->sampleTime;
449      simnfo->thermalTime = simnfo->sampleTime;
450    }
# Line 326 | Line 460 | void SimSetup::createSim( void ){
460    // check for the temperature set flag
461  
462    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
463 <  
464 <    
463 >
464 >
465    // make the longe range forces and the integrator
466 <  
466 >
467    new AllLong( simnfo );
468 <      
468 >
469    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
470    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
471    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
472   }
473  
474   void SimSetup::makeAtoms( void ){
475 <  
475 >
476    int i, j, k, index;
477    double ux, uy, uz, uSqr, u;
478    AtomStamp* current_atom;
479    DirectionalAtom* dAtom;
480    int molIndex, molStart, molEnd, nMemb;
347 –  
481  
482 +
483    molIndex = 0;
484    index = 0;
485    for( i=0; i<n_components; i++ ){
486 <    
486 >
487      for( j=0; j<components_nmol[i]; j++ ){
488 <      
488 >
489        molStart = index;
490        nMemb = comp_stamps[i]->getNAtoms();
491        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
492 <        
492 >
493          current_atom = comp_stamps[i]->getAtom( k );
494 <        if( current_atom->haveOrientation() ){
494 >        if( current_atom->haveOrientation() ){
495  
496 <          dAtom = new DirectionalAtom;
496 >          dAtom = new DirectionalAtom(index);
497            simnfo->n_oriented++;
498            the_atoms[index] = dAtom;
499 <      
499 >
500            ux = current_atom->getOrntX();
501            uy = current_atom->getOrntY();
502            uz = current_atom->getOrntZ();
503 <          
503 >
504            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
505 <          
505 >
506            u = sqrt( uSqr );
507            ux = ux / u;
508            uy = uy / u;
509            uz = uz / u;
510 <          
510 >
511            dAtom->setSUx( ux );
512            dAtom->setSUy( uy );
513            dAtom->setSUz( uz );
514          }
515          else{
516 <          the_atoms[index] = new GeneralAtom;
516 >          the_atoms[index] = new GeneralAtom(index);
517          }
518          the_atoms[index]->setType( current_atom->getType() );
519          the_atoms[index]->setIndex( index );
520 <        
520 >
521          // increment the index and repeat;
522          index++;
523        }
524 <      
524 >
525        molEnd = index -1;
526        the_molecules[molIndex].setNMembers( nMemb );
527        the_molecules[molIndex].setStartAtom( molStart );
# Line 396 | Line 530 | void SimSetup::makeAtoms( void ){
530  
531      }
532    }
533 <  
533 >
534    the_ff->initializeAtoms();
535   }
536  
# Line 410 | Line 544 | void SimSetup::makeBonds( void ){
544    index = 0;
545    offset = 0;
546    for( i=0; i<n_components; i++ ){
547 <    
547 >
548      for( j=0; j<components_nmol[i]; j++ ){
549 <      
549 >
550        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
551 <        
551 >
552          current_bond = comp_stamps[i]->getBond( k );
553          the_bonds[index].a = current_bond->getA() + offset;
554          the_bonds[index].b = current_bond->getB() + offset;
# Line 428 | Line 562 | void SimSetup::makeBonds( void ){
562        offset += comp_stamps[i]->getNAtoms();
563      }
564    }
565 <  
565 >
566    the_ff->initializeBonds( the_bonds );
567   }
568  
# Line 442 | Line 576 | void SimSetup::makeBends( void ){
576    index = 0;
577    offset = 0;
578    for( i=0; i<n_components; i++ ){
579 <    
579 >
580      for( j=0; j<components_nmol[i]; j++ ){
581 <      
581 >
582        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
583 <        
583 >
584          current_bend = comp_stamps[i]->getBend( k );
585          the_bends[index].a = current_bend->getA() + offset;
586          the_bends[index].b = current_bend->getB() + offset;
# Line 461 | Line 595 | void SimSetup::makeBends( void ){
595        offset += comp_stamps[i]->getNAtoms();
596      }
597    }
598 <  
598 >
599    the_ff->initializeBends( the_bends );
600   }
601  
# Line 475 | Line 609 | void SimSetup::makeTorsions( void ){
609    index = 0;
610    offset = 0;
611    for( i=0; i<n_components; i++ ){
612 <    
612 >
613      for( j=0; j<components_nmol[i]; j++ ){
614 <      
614 >
615        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
616 <        
616 >
617          current_torsion = comp_stamps[i]->getTorsion( k );
618          the_torsions[index].a = current_torsion->getA() + offset;
619          the_torsions[index].b = current_torsion->getB() + offset;
# Line 495 | Line 629 | void SimSetup::makeTorsions( void ){
629        offset += comp_stamps[i]->getNAtoms();
630      }
631    }
632 <  
632 >
633    the_ff->initializeTorsions( the_torsions );
634   }
635  
# Line 526 | Line 660 | void SimSetup::initFromBass( void ){
660      n_per_extra = (int)ceil( temp1 );
661  
662      if( n_per_extra > 4){
663 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
664 <      exit(8);
663 >      sprintf( painCave.errMsg,
664 >               "SimSetup error. There has been an error in constructing"
665 >               " the non-complete lattice.\n" );
666 >      painCave.isFatal = 1;
667 >      simError();
668      }
669    }
670    else{
# Line 536 | Line 673 | void SimSetup::initFromBass( void ){
673      celly = simnfo->box_y / temp3;
674      cellz = simnfo->box_z / temp3;
675    }
676 <  
676 >
677    current_mol = 0;
678    current_comp_mol = 0;
679    current_comp = 0;
680    current_atom_ndx = 0;
681 <  
681 >
682    for( i=0; i < n_cells ; i++ ){
683      for( j=0; j < n_cells; j++ ){
684        for( k=0; k < n_cells; k++ ){
685 <        
685 >
686          makeElement( i * cellx,
687                       j * celly,
688                       k * cellz );
689 <        
689 >
690          makeElement( i * cellx + 0.5 * cellx,
691                       j * celly + 0.5 * celly,
692                       k * cellz );
693 <        
693 >
694          makeElement( i * cellx,
695                       j * celly + 0.5 * celly,
696                       k * cellz + 0.5 * cellz );
697 <        
697 >
698          makeElement( i * cellx + 0.5 * cellx,
699                       j * celly,
700                       k * cellz + 0.5 * cellz );
# Line 567 | Line 704 | void SimSetup::initFromBass( void ){
704  
705    if( have_extra ){
706      done = 0;
707 <    
707 >
708      int start_ndx;
709      for( i=0; i < (n_cells+1) && !done; i++ ){
710        for( j=0; j < (n_cells+1) && !done; j++ ){
711 <        
711 >
712          if( i < n_cells ){
713 <          
713 >
714            if( j < n_cells ){
715              start_ndx = n_cells;
716            }
717            else start_ndx = 0;
718          }
719          else start_ndx = 0;
720 <        
720 >
721          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
722 <          
722 >
723            makeElement( i * cellx,
724                         j * celly,
725                         k * cellz );
726            done = ( current_mol >= tot_nmol );
727 <          
727 >
728            if( !done && n_per_extra > 1 ){
729              makeElement( i * cellx + 0.5 * cellx,
730                           j * celly + 0.5 * celly,
731                           k * cellz );
732              done = ( current_mol >= tot_nmol );
733            }
734 <          
734 >
735            if( !done && n_per_extra > 2){
736              makeElement( i * cellx,
737                           j * celly + 0.5 * celly,
738                           k * cellz + 0.5 * cellz );
739              done = ( current_mol >= tot_nmol );
740            }
741 <          
741 >
742            if( !done && n_per_extra > 3){
743              makeElement( i * cellx + 0.5 * cellx,
744                           j * celly,
# Line 612 | Line 749 | void SimSetup::initFromBass( void ){
749        }
750      }
751    }
752 <  
753 <  
752 >
753 >
754    for( i=0; i<simnfo->n_atoms; i++ ){
755      simnfo->atoms[i]->set_vx( 0.0 );
756      simnfo->atoms[i]->set_vy( 0.0 );
# Line 629 | Line 766 | void SimSetup::makeElement( double x, double y, double
766    double rotMat[3][3];
767  
768    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
769 <    
769 >
770      current_atom = comp_stamps[current_comp]->getAtom( k );
771      if( !current_atom->havePosition() ){
772 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
773 <                << ", atom " << current_atom->getType()
774 <                << " does not have a position specified.\n"
775 <                << "The initialization routine is unable to give a start"
776 <                << " position.\n";
777 <      exit(8);
772 >      sprintf( painCave.errMsg,
773 >               "SimSetup:initFromBass error.\n"
774 >               "\tComponent %s, atom %s does not have a position specified.\n"
775 >               "\tThe initialization routine is unable to give a start"
776 >               " position.\n",
777 >               comp_stamps[current_comp]->getID(),
778 >               current_atom->getType() );
779 >      painCave.isFatal = 1;
780 >      simError();
781      }
782 <    
782 >
783      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
784      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
785      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
786 <    
786 >
787      if( the_atoms[current_atom_ndx]->isDirectional() ){
788 <      
788 >
789        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
790 <      
790 >
791        rotMat[0][0] = 1.0;
792        rotMat[0][1] = 0.0;
793        rotMat[0][2] = 0.0;
# Line 665 | Line 805 | void SimSetup::makeElement( double x, double y, double
805  
806      current_atom_ndx++;
807    }
808 <  
808 >
809    current_mol++;
810    current_comp_mol++;
811  
812    if( current_comp_mol >= components_nmol[current_comp] ){
813 <    
813 >
814      current_comp_mol = 0;
815      current_comp++;
816    }

Diff Legend

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