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 11 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
Revision 202 by mmeineke, Tue Dec 10 21:41:26 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;
75    Globals* the_globals;
76 <  int i;
76 >  int i, j;
77  
78    // get the stamps and globals;
79    the_stamps = 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 +  // also extract the used stamps out into a separate linked list
165  
166 +  simnfo->nComponents = n_components;
167 +  simnfo->componentsNmol = components_nmol;
168 +  simnfo->compStamps = comp_stamps;
169 +  simnfo->headStamp = new LinkedMolStamp();
170 +  
171 +  char* id;
172 +  LinkedMolStamp* headStamp = simnfo->headStamp;
173 +  LinkedMolStamp* currentStamp = NULL;
174    for( i=0; i<n_components; i++ ){
175  
176 <    comp_stamps[i] =
177 <      the_stamps->getMolecule( the_components[i]->getType() );
176 >    id = the_components[i]->getType();
177 >    comp_stamps[i] = NULL;
178 >    
179 >    // check to make sure the component isn't already in the list
180 >
181 >    comp_stamps[i] = headStamp->match( id );
182 >    if( comp_stamps[i] == NULL ){
183 >      
184 >      // extract the component from the list;
185 >      
186 >      currentStamp = the_stamps->extractMolStamp( id );
187 >      if( currentStamp == NULL ){
188 >        sprintf( painCave.errMsg,
189 >                 "SimSetup error: Component \"%s\" was not found in the "
190 >                 "list of declared molecules\n"
191 >                 id );
192 >        painCave.isFatal = 1;
193 >        simError();
194 >      }
195 >      
196 >      headStamp->add( currentStamp );
197 >      comp_stamps[i] = headStamp->match( id );
198 >    }
199    }
200  
201 + #ifdef IS_MPI
202 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
203 +  MPIcheckPoint();
204 + #endif // is_mpi
205    
206  
207 +
208 +
209    // caclulate the number of atoms, bonds, bends and torsions
210  
211    tot_atoms = 0;
# Line 111 | Line 214 | void SimSetup::createSim( void ){
214    tot_torsions = 0;
215    for( i=0; i<n_components; i++ ){
216      
217 <    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
218 <    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
219 <    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
217 >    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
218 >    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
219 >    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
220      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
221    }
222 <  
222 >
223    tot_SRI = tot_bonds + tot_bends + tot_torsions;
224 <  
224 >
225    simnfo->n_atoms = tot_atoms;
226    simnfo->n_bonds = tot_bonds;
227    simnfo->n_bends = tot_bends;
228    simnfo->n_torsions = tot_torsions;
229    simnfo->n_SRI = tot_SRI;
230 +  simnfo->n_mol = tot_nmol;
231  
128 –  // create the atom and short range interaction arrays
232    
233 <  the_atoms = new Atom*[tot_atoms];
234 <  //  the_molecules = new Molecule[tot_nmol];
233 > #ifdef IS_MPI
234 >
235 >  // divide the molecules among processors here.
236    
237 +  mpiSimulation* mpiSim = new mpiSimulation( simnfo );
238    
239 <  if( tot_SRI ){
240 <    the_sris = new SRI*[tot_SRI];
241 <    the_excludes = new ex_pair[tot_SRI];
239 >  mpiSim->divideLabor();
240 >
241 >  // set up the local variables
242 >  
243 >  int localMol;
244 >  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
245 >  
246 >  localMol = 0;
247 >  local_atoms = 0;
248 >  local_bonds = 0;
249 >  local_bends = 0;
250 >  local_torsions = 0;
251 >  for( i=0; i<n_components; i++ ){
252 >
253 >    for( j=0; j<components_nmol[i]; j++ ){
254 >      
255 >      if( mpiSim->getMyMolStart() <= j &&
256 >          j <= mpiSim->getMyMolEnd() ){
257 >        
258 >        local_atoms +=    comp_stamps[i]->getNAtoms();
259 >        local_bonds +=    comp_stamps[i]->getNBonds();
260 >        local_bends +=    comp_stamps[i]->getNBends();
261 >        local_torsions += comp_stamps[i]->getNTorsions();
262 >        localMol++;
263 >      }      
264 >    }
265    }
266  
267 +  
268 +
269 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
270 +  
271 +
272 + #endif // is_mpi
273 +  
274 +
275 +  // create the atom and short range interaction arrays
276 +
277 +  Atom::createArrays(simnfo->n_atoms);
278 +  the_atoms = new Atom*[simnfo->n_atoms];
279 +  the_molecules = new Molecule[simnfo->n_mol];
280 +
281 +
282 +  if( simnfo->n_SRI ){
283 +    the_sris = new SRI*[simnfo->n_SRI];
284 +    the_excludes = new ex_pair[simnfo->n_SRI];
285 +  }
286 +
287    // set the arrays into the SimInfo object
288  
289    simnfo->atoms = the_atoms;
# Line 143 | Line 291 | void SimSetup::createSim( void ){
291    simnfo->n_exclude = tot_SRI;
292    simnfo->excludes = the_excludes;
293  
294 +
295    // initialize the arrays
296 <  
296 >
297    the_ff->setSimInfo( simnfo );
298 <    
298 >
299    makeAtoms();
300  
301    if( tot_bonds ){
# Line 161 | Line 310 | void SimSetup::createSim( void ){
310      makeTorsions();
311    }
312  
164 –  //  makeMolecules();
313  
314    // get some of the tricky things that may still be in the globals
315  
316    if( simnfo->n_dipoles ){
317  
318      if( !the_globals->haveRRF() ){
319 <      std::cerr << "SimSetup Error, system has dipoles, but no rRF was set.\n";
320 <      exit(8);
319 >      sprintf( painCave.errMsg,
320 >               "SimSetup Error, system has dipoles, but no rRF was set.\n");
321 >      painCave.isFatal = 1;
322 >      simError();
323      }
324      if( !the_globals->haveDielectric() ){
325 <      std::cerr << "SimSetup Error, system has dipoles, but no"
326 <                << " dielectric was set.\n";
327 <      exit(8);
325 >      sprintf( painCave.errMsg,
326 >               "SimSetup Error, system has dipoles, but no"
327 >               " dielectric was set.\n" );
328 >      painCave.isFatal = 1;
329 >      simError();
330      }
331  
332      simnfo->rRF        = the_globals->getRRF();
333      simnfo->dielectric = the_globals->getDielectric();
334    }
335  
336 + #ifdef IS_MPI
337 +  strcpy( checkPointMsg, "rRf and dielectric check out" );
338 +  MPIcheckPoint();
339 + #endif // is_mpi
340 +  
341    if( the_globals->haveBox() ){
342      simnfo->box_x = the_globals->getBox();
343      simnfo->box_y = the_globals->getBox();
344      simnfo->box_z = the_globals->getBox();
345    }
346    else if( the_globals->haveDensity() ){
347 <    
347 >
348      double vol;
349      vol = (double)tot_nmol / the_globals->getDensity();
350      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 196 | Line 353 | void SimSetup::createSim( void ){
353    }
354    else{
355      if( !the_globals->haveBoxX() ){
356 <      std::cerr << "SimSetup error, no periodic BoxX size given.\n";
357 <      exit(8);
356 >      sprintf( painCave.errMsg,
357 >               "SimSetup error, no periodic BoxX size given.\n" );
358 >      painCave.isFatal = 1;
359 >      simError();
360      }
361      simnfo->box_x = the_globals->getBoxX();
362  
363      if( !the_globals->haveBoxY() ){
364 <      std::cerr << "SimSetup error, no periodic BoxY size given.\n";
365 <      exit(8);
364 >      sprintf( painCave.errMsg,
365 >               "SimSetup error, no periodic BoxY size given.\n" );
366 >      painCave.isFatal = 1;
367 >      simError();
368      }
369      simnfo->box_y = the_globals->getBoxY();
370  
371      if( !the_globals->haveBoxZ() ){
372 <      std::cerr << "SimSetup error, no periodic BoxZ size given.\n";
373 <      exit(8);
372 >      sprintf( painCave.errMsg,
373 >               "SimSetup error, no periodic BoxZ size given.\n" );
374 >      painCave.isFatal = 1;
375 >      simError();
376      }
377      simnfo->box_z = the_globals->getBoxZ();
378    }
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
379  
380 <    delete fileInit;    
381 <  }
382 <  else{
383 <    initFromBass();
227 <  }
380 > #ifdef IS_MPI
381 >  strcpy( checkPointMsg, "Box size set up" );
382 >  MPIcheckPoint();
383 > #endif // is_mpi
384  
385 <  if( the_globals->haveFinalConfig() ){
386 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
387 <  }
388 <  else{
389 <    strcpy( simnfo->finalName, inFileName );
390 <    char* endTest;
391 <    int nameLength = strlen( simnfo->finalName );
392 <    endTest = &(simnfo->finalName[nameLength - 5]);
393 <    if( !strcmp( endTest, ".bass" ) ){
394 <      strcpy( endTest, ".eor" );
395 <    }
396 <    else if( !strcmp( endTest, ".BASS" ) ){
397 <      strcpy( endTest, ".eor" );
385 >
386 >
387 > //   if( the_globals->haveInitialConfig() ){
388 > //        InitializeFromFile* fileInit;
389 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
390 >
391 > //     fileInit->read_xyz( simnfo ); // default velocities on
392 >
393 > //     delete fileInit;
394 > //   }
395 > //   else{
396 >
397 > #ifdef IS_MPI
398 >
399 >  // no init from bass
400 >  
401 >  sprintf( painCave.errMsg,
402 >           "Cannot intialize a parallel simulation without an initial configuration file.\n" );
403 >  painCave.isFatal;
404 >  simError();
405 >  
406 > #else
407 >
408 >  initFromBass();
409 >
410 > #endif // is_mpi
411 >
412 > #ifdef IS_MPI
413 >  strcpy( checkPointMsg, "Successfully read in the initial configuration" );
414 >  MPIcheckPoint();
415 > #endif // is_mpi
416 >
417 >
418 >  
419 >
420 >  
421 >  //   }
422 >  
423 > #ifdef IS_MPI
424 >  if( worldRank == 0 ){
425 > #endif // is_mpi
426 >    
427 >    if( the_globals->haveFinalConfig() ){
428 >      strcpy( simnfo->finalName, the_globals->getFinalConfig() );
429      }
430      else{
431 <      endTest = &(simnfo->finalName[nameLength - 4]);
432 <      if( !strcmp( endTest, ".bss" ) ){
431 >      strcpy( simnfo->finalName, inFileName );
432 >      char* endTest;
433 >      int nameLength = strlen( simnfo->finalName );
434 >      endTest = &(simnfo->finalName[nameLength - 5]);
435 >      if( !strcmp( endTest, ".bass" ) ){
436          strcpy( endTest, ".eor" );
437        }
438 <      else if( !strcmp( endTest, ".mdl" ) ){
438 >      else if( !strcmp( endTest, ".BASS" ) ){
439          strcpy( endTest, ".eor" );
440        }
441        else{
442 <        strcat( simnfo->finalName, ".eor" );
442 >        endTest = &(simnfo->finalName[nameLength - 4]);
443 >        if( !strcmp( endTest, ".bss" ) ){
444 >          strcpy( endTest, ".eor" );
445 >        }
446 >        else if( !strcmp( endTest, ".mdl" ) ){
447 >          strcpy( endTest, ".eor" );
448 >        }
449 >        else{
450 >          strcat( simnfo->finalName, ".eor" );
451 >        }
452        }
453      }
255 –  }
454      
455 <  // make the sample and status out names
456 <
457 <  strcpy( simnfo->sampleName, inFileName );
458 <  char* endTest;
459 <  int nameLength = strlen( simnfo->sampleName );
460 <  endTest = &(simnfo->sampleName[nameLength - 5]);
461 <  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" ) ){
455 >    // make the sample and status out names
456 >    
457 >    strcpy( simnfo->sampleName, inFileName );
458 >    char* endTest;
459 >    int nameLength = strlen( simnfo->sampleName );
460 >    endTest = &(simnfo->sampleName[nameLength - 5]);
461 >    if( !strcmp( endTest, ".bass" ) ){
462        strcpy( endTest, ".dump" );
463      }
464 <    else if( !strcmp( endTest, ".mdl" ) ){
464 >    else if( !strcmp( endTest, ".BASS" ) ){
465        strcpy( endTest, ".dump" );
466      }
467      else{
468 <      strcat( simnfo->sampleName, ".dump" );
468 >      endTest = &(simnfo->sampleName[nameLength - 4]);
469 >      if( !strcmp( endTest, ".bss" ) ){
470 >        strcpy( endTest, ".dump" );
471 >      }
472 >      else if( !strcmp( endTest, ".mdl" ) ){
473 >        strcpy( endTest, ".dump" );
474 >      }
475 >      else{
476 >        strcat( simnfo->sampleName, ".dump" );
477 >      }
478      }
479 <  }
480 <  
481 <  strcpy( simnfo->statusName, inFileName );
482 <  nameLength = strlen( simnfo->statusName );
483 <  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" ) ){
479 >    
480 >    strcpy( simnfo->statusName, inFileName );
481 >    nameLength = strlen( simnfo->statusName );
482 >    endTest = &(simnfo->statusName[nameLength - 5]);
483 >    if( !strcmp( endTest, ".bass" ) ){
484        strcpy( endTest, ".stat" );
485      }
486 <    else if( !strcmp( endTest, ".mdl" ) ){
486 >    else if( !strcmp( endTest, ".BASS" ) ){
487        strcpy( endTest, ".stat" );
488      }
489      else{
490 <      strcat( simnfo->statusName, ".stat" );
490 >      endTest = &(simnfo->statusName[nameLength - 4]);
491 >      if( !strcmp( endTest, ".bss" ) ){
492 >        strcpy( endTest, ".stat" );
493 >      }
494 >      else if( !strcmp( endTest, ".mdl" ) ){
495 >        strcpy( endTest, ".stat" );
496 >      }
497 >      else{
498 >        strcat( simnfo->statusName, ".stat" );
499 >      }
500      }
501 +    
502 + #ifdef IS_MPI
503    }
504 + #endif // is_mpi
505    
506    // set the status, sample, and themal kick times
507 <
507 >  
508    if( the_globals->haveSampleTime() ){
509 <    simnfo->sampleTime = the_globals->getSampleTime();
509 >    simnfo->sampleTime = the_globals->getSampleTime();
510      simnfo->statusTime = simnfo->sampleTime;
511      simnfo->thermalTime = simnfo->sampleTime;
512    }
513    else{
514 <    simnfo->sampleTime = the_globals->getRunTime();
514 >    simnfo->sampleTime = the_globals->getRunTime();
515      simnfo->statusTime = simnfo->sampleTime;
516      simnfo->thermalTime = simnfo->sampleTime;
517    }
# Line 325 | Line 527 | void SimSetup::createSim( void ){
527    // check for the temperature set flag
528  
529    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
530 <  
531 <    
530 >
531 >
532    // make the longe range forces and the integrator
533 <  
533 >
534    new AllLong( simnfo );
535 <      
535 >
536    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
537    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
538    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
539   }
540  
541   void SimSetup::makeAtoms( void ){
542 <  
542 >
543    int i, j, k, index;
544    double ux, uy, uz, uSqr, u;
545    AtomStamp* current_atom;
546    DirectionalAtom* dAtom;
547 +  int molIndex, molStart, molEnd, nMemb;
548  
549 +
550 +  molIndex = 0;
551    index = 0;
552    for( i=0; i<n_components; i++ ){
553 <    
553 >
554      for( j=0; j<components_nmol[i]; j++ ){
555 <      
555 >
556 >      molStart = index;
557 >      nMemb = comp_stamps[i]->getNAtoms();
558        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
559 <        
559 >
560          current_atom = comp_stamps[i]->getAtom( k );
561 <        if( current_atom->haveOrientation() ){
561 >        if( current_atom->haveOrientation() ){
562  
563 <          dAtom = new DirectionalAtom;
563 >          dAtom = new DirectionalAtom(index);
564            simnfo->n_oriented++;
565            the_atoms[index] = dAtom;
566 <      
566 >
567            ux = current_atom->getOrntX();
568            uy = current_atom->getOrntY();
569            uz = current_atom->getOrntZ();
570 <          
570 >
571            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
572 <          
572 >
573            u = sqrt( uSqr );
574            ux = ux / u;
575            uy = uy / u;
576            uz = uz / u;
577 <          
577 >
578            dAtom->setSUx( ux );
579            dAtom->setSUy( uy );
580            dAtom->setSUz( uz );
581          }
582          else{
583 <          the_atoms[index] = new GeneralAtom;
583 >          the_atoms[index] = new GeneralAtom(index);
584          }
585          the_atoms[index]->setType( current_atom->getType() );
586          the_atoms[index]->setIndex( index );
587 <        
587 >
588          // increment the index and repeat;
589          index++;
590        }
591 +
592 +      molEnd = index -1;
593 +      the_molecules[molIndex].setNMembers( nMemb );
594 +      the_molecules[molIndex].setStartAtom( molStart );
595 +      the_molecules[molIndex].setEndAtom( molEnd );
596 +      the_molecules[molIndex].setStampID( i );
597 +      molIndex++;
598 +
599      }
600    }
601 <  
601 >
602    the_ff->initializeAtoms();
603   }
604  
# Line 397 | Line 612 | void SimSetup::makeBonds( void ){
612    index = 0;
613    offset = 0;
614    for( i=0; i<n_components; i++ ){
615 <    
615 >
616      for( j=0; j<components_nmol[i]; j++ ){
617 <      
617 >
618        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
619 <        
619 >
620          current_bond = comp_stamps[i]->getBond( k );
621          the_bonds[index].a = current_bond->getA() + offset;
622          the_bonds[index].b = current_bond->getB() + offset;
# Line 415 | Line 630 | void SimSetup::makeBonds( void ){
630        offset += comp_stamps[i]->getNAtoms();
631      }
632    }
633 <  
633 >
634    the_ff->initializeBonds( the_bonds );
635   }
636  
# Line 429 | Line 644 | void SimSetup::makeBends( void ){
644    index = 0;
645    offset = 0;
646    for( i=0; i<n_components; i++ ){
647 <    
647 >
648      for( j=0; j<components_nmol[i]; j++ ){
649 <      
649 >
650        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
651 <        
651 >
652          current_bend = comp_stamps[i]->getBend( k );
653          the_bends[index].a = current_bend->getA() + offset;
654          the_bends[index].b = current_bend->getB() + offset;
# Line 448 | Line 663 | void SimSetup::makeBends( void ){
663        offset += comp_stamps[i]->getNAtoms();
664      }
665    }
666 <  
666 >
667    the_ff->initializeBends( the_bends );
668   }
669  
# Line 462 | Line 677 | void SimSetup::makeTorsions( void ){
677    index = 0;
678    offset = 0;
679    for( i=0; i<n_components; i++ ){
680 <    
680 >
681      for( j=0; j<components_nmol[i]; j++ ){
682 <      
682 >
683        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
684 <        
684 >
685          current_torsion = comp_stamps[i]->getTorsion( k );
686          the_torsions[index].a = current_torsion->getA() + offset;
687          the_torsions[index].b = current_torsion->getB() + offset;
# Line 482 | Line 697 | void SimSetup::makeTorsions( void ){
697        offset += comp_stamps[i]->getNAtoms();
698      }
699    }
700 <  
700 >
701    the_ff->initializeTorsions( the_torsions );
702   }
703  
489 – void SimSetup::makeMolecules( void ){
490 –
491 –  //empy for now
492 – }
493 –
704   void SimSetup::initFromBass( void ){
705  
706    int i, j, k;
# Line 518 | Line 728 | void SimSetup::initFromBass( void ){
728      n_per_extra = (int)ceil( temp1 );
729  
730      if( n_per_extra > 4){
731 <      std::cerr << "THere has been an error in constructing the non-complete lattice.\n";
732 <      exit(8);
731 >      sprintf( painCave.errMsg,
732 >               "SimSetup error. There has been an error in constructing"
733 >               " the non-complete lattice.\n" );
734 >      painCave.isFatal = 1;
735 >      simError();
736      }
737    }
738    else{
# Line 528 | Line 741 | void SimSetup::initFromBass( void ){
741      celly = simnfo->box_y / temp3;
742      cellz = simnfo->box_z / temp3;
743    }
744 <  
744 >
745    current_mol = 0;
746    current_comp_mol = 0;
747    current_comp = 0;
748    current_atom_ndx = 0;
749 <  
749 >
750    for( i=0; i < n_cells ; i++ ){
751      for( j=0; j < n_cells; j++ ){
752        for( k=0; k < n_cells; k++ ){
753 <        
753 >
754          makeElement( i * cellx,
755                       j * celly,
756                       k * cellz );
757 <        
757 >
758          makeElement( i * cellx + 0.5 * cellx,
759                       j * celly + 0.5 * celly,
760                       k * cellz );
761 <        
761 >
762          makeElement( i * cellx,
763                       j * celly + 0.5 * celly,
764                       k * cellz + 0.5 * cellz );
765 <        
765 >
766          makeElement( i * cellx + 0.5 * cellx,
767                       j * celly,
768                       k * cellz + 0.5 * cellz );
# Line 559 | Line 772 | void SimSetup::initFromBass( void ){
772  
773    if( have_extra ){
774      done = 0;
775 <    
775 >
776      int start_ndx;
777      for( i=0; i < (n_cells+1) && !done; i++ ){
778        for( j=0; j < (n_cells+1) && !done; j++ ){
779 <        
779 >
780          if( i < n_cells ){
781 <          
781 >
782            if( j < n_cells ){
783              start_ndx = n_cells;
784            }
785            else start_ndx = 0;
786          }
787          else start_ndx = 0;
788 <        
788 >
789          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
790 <          
790 >
791            makeElement( i * cellx,
792                         j * celly,
793                         k * cellz );
794            done = ( current_mol >= tot_nmol );
795 <          
795 >
796            if( !done && n_per_extra > 1 ){
797              makeElement( i * cellx + 0.5 * cellx,
798                           j * celly + 0.5 * celly,
799                           k * cellz );
800              done = ( current_mol >= tot_nmol );
801            }
802 <          
802 >
803            if( !done && n_per_extra > 2){
804              makeElement( i * cellx,
805                           j * celly + 0.5 * celly,
806                           k * cellz + 0.5 * cellz );
807              done = ( current_mol >= tot_nmol );
808            }
809 <          
809 >
810            if( !done && n_per_extra > 3){
811              makeElement( i * cellx + 0.5 * cellx,
812                           j * celly,
# Line 604 | Line 817 | void SimSetup::initFromBass( void ){
817        }
818      }
819    }
820 <  
821 <  
820 >
821 >
822    for( i=0; i<simnfo->n_atoms; i++ ){
823      simnfo->atoms[i]->set_vx( 0.0 );
824      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 834 | void SimSetup::makeElement( double x, double y, double
834    double rotMat[3][3];
835  
836    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
837 <    
837 >
838      current_atom = comp_stamps[current_comp]->getAtom( k );
839      if( !current_atom->havePosition() ){
840 <      std::cerr << "Component " << comp_stamps[current_comp]->getID()
841 <                << ", atom " << current_atom->getType()
842 <                << " does not have a position specified.\n"
843 <                << "The initialization routine is unable to give a start"
844 <                << " position.\n";
845 <      exit(8);
840 >      sprintf( painCave.errMsg,
841 >               "SimSetup:initFromBass error.\n"
842 >               "\tComponent %s, atom %s does not have a position specified.\n"
843 >               "\tThe initialization routine is unable to give a start"
844 >               " position.\n",
845 >               comp_stamps[current_comp]->getID(),
846 >               current_atom->getType() );
847 >      painCave.isFatal = 1;
848 >      simError();
849      }
850 <    
850 >
851      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
852      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
853      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
854 <    
854 >
855      if( the_atoms[current_atom_ndx]->isDirectional() ){
856 <      
856 >
857        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
858 <      
858 >
859        rotMat[0][0] = 1.0;
860        rotMat[0][1] = 0.0;
861        rotMat[0][2] = 0.0;
# Line 657 | Line 873 | void SimSetup::makeElement( double x, double y, double
873  
874      current_atom_ndx++;
875    }
876 <  
876 >
877    current_mol++;
878    current_comp_mol++;
879  
880    if( current_comp_mol >= components_nmol[current_comp] ){
881 <    
881 >
882      current_comp_mol = 0;
883      current_comp++;
884    }

Diff Legend

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