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 164 by mmeineke, Tue Nov 5 22:04:46 2002 UTC

# Line 7 | Line 7
7   #include "LRI.hpp"
8   #include "Integrator.hpp"
9  
10 + #ifdef IS_MPI
11 + #include "mpiBASS.h"
12 + #include "bassDiag.hpp"
13 + #endif
14 +
15   SimSetup::SimSetup(){
16    stamps = new MakeStamps();
17    globals = new Globals();
18 +  
19 + #ifdef IS_MPI
20 +  strcpy( checkPointMsg, "SimSetup creation successful" );
21 +  MPIcheckPoint();
22 + #endif IS_MPI
23   }
24  
25   SimSetup::~SimSetup(){
# Line 19 | Line 29 | void SimSetup::parseFile( char* fileName ){
29  
30   void SimSetup::parseFile( char* fileName ){
31  
32 <  inFileName = fileName;
33 <  set_interface_stamps( stamps, globals );
34 <  yacc_BASS( fileName );
32 > #ifdef IS_MPI
33 >  if( worldRank == 0 ){
34 > #endif // is_mpi
35 >    
36 >    inFileName = fileName;
37 >    set_interface_stamps( stamps, globals );
38 >    
39 > #ifdef IS_MPI
40 >    mpiEventInit();
41 > #endif
42 >
43 >    yacc_BASS( fileName );
44 >
45 > #ifdef IS_MPI
46 >    throwMPIEvent(NULL);
47 >  }
48 >  else recieveParse();
49 > #endif
50 >
51   }
52  
53 + #ifdef IS_MPI
54 + void SimSetup::receiveParse(void){
55 +
56 +    set_interface_stamps( stamps, globals );
57 +    mpiEventInit();
58 +    MPIcheckPoint();
59 +    mpiEventLoop();
60 +
61 + }
62 +
63 +
64 + void SimSetup::testMe(void){
65 +  bassDiag* dumpMe = new bassDiag(globals,stamps);
66 +  dumpMe->dumpStamps();
67 +  delete dumpMe;
68 + }
69 + #endif
70   void SimSetup::createSim( void ){
71  
72    MakeStamps *the_stamps;
# Line 43 | Line 86 | void SimSetup::createSim( void ){
86    n_components = the_globals->getNComponents();
87    strcpy( force_field, the_globals->getForceField() );
88    strcpy( ensemble, the_globals->getEnsemble() );
89 <  
89 >
90    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
91    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
92    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
93    else{
94 <    std::cerr<< "SimSetup Error. Unrecognized force field -> "
94 >    std::cerr<< "SimSetup Error. Unrecognized force field -> "
95               << force_field << "\n";
96      exit(8);
97    }
# Line 57 | Line 100 | void SimSetup::createSim( void ){
100    the_components = the_globals->getComponents();
101    components_nmol = new int[n_components];
102    comp_stamps = new MoleculeStamp*[n_components];
103 <  
103 >
104    if( !the_globals->haveNMol() ){
105 <    // we don't have the total number of molecules, so we assume it is
105 >    // we don't have the total number of molecules, so we assume it is
106      // given in each component
107  
108      tot_nmol = 0;
109      for( i=0; i<n_components; i++ ){
110 <      
110 >
111        if( !the_components[i]->haveNMol() ){
112          // we have a problem
113          std::cerr << "SimSetup Error. No global NMol or component NMol"
# Line 78 | Line 121 | void SimSetup::createSim( void ){
121    }
122    else{
123      std::cerr << "NOT A SUPPORTED FEATURE\n";
124 <    
124 >
125   //     tot_nmol = the_globals->getNMol();
126 <    
126 >
127   //     //we have the total number of molecules, now we check for molfractions
128   //     for( i=0; i<n_components; i++ ){
129 <      
129 >
130   //       if( !the_components[i]->haveMolFraction() ){
131 <        
131 >
132   //      if( !the_components[i]->haveNMol() ){
133   //        //we have a problem
134   //        std::cerr << "SimSetup error. Neither molFraction nor "
# Line 97 | Line 140 | void SimSetup::createSim( void ){
140  
141    for( i=0; i<n_components; i++ ){
142  
143 <    comp_stamps[i] =
143 >    comp_stamps[i] =
144        the_stamps->getMolecule( the_components[i]->getType() );
145    }
146  
104 –  
147  
148 +
149    // caclulate the number of atoms, bonds, bends and torsions
150  
151    tot_atoms = 0;
# Line 110 | Line 153 | void SimSetup::createSim( void ){
153    tot_bends = 0;
154    tot_torsions = 0;
155    for( i=0; i<n_components; i++ ){
156 <    
156 >
157      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
158      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
159      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
160      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
161    }
162 <  
162 >
163    tot_SRI = tot_bonds + tot_bends + tot_torsions;
164 <  
164 >
165    simnfo->n_atoms = tot_atoms;
166    simnfo->n_bonds = tot_bonds;
167    simnfo->n_bends = tot_bends;
# Line 126 | Line 169 | void SimSetup::createSim( void ){
169    simnfo->n_SRI = tot_SRI;
170  
171    // create the atom and short range interaction arrays
172 <  
172 >
173    the_atoms = new Atom*[tot_atoms];
174 <  //  the_molecules = new Molecule[tot_nmol];
175 <  
176 <  
174 >  Atom::createArrays(tot_atoms);
175 >  the_molecules = new Molecule[tot_nmol];
176 >
177 >
178    if( tot_SRI ){
179      the_sris = new SRI*[tot_SRI];
180      the_excludes = new ex_pair[tot_SRI];
# Line 143 | Line 187 | void SimSetup::createSim( void ){
187    simnfo->n_exclude = tot_SRI;
188    simnfo->excludes = the_excludes;
189  
190 +
191    // initialize the arrays
192 <  
192 >
193    the_ff->setSimInfo( simnfo );
194 <    
194 >
195    makeAtoms();
196  
197    if( tot_bonds ){
# Line 187 | Line 232 | void SimSetup::createSim( void ){
232      simnfo->box_z = the_globals->getBox();
233    }
234    else if( the_globals->haveDensity() ){
235 <    
235 >
236      double vol;
237      vol = (double)tot_nmol / the_globals->getDensity();
238      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 213 | Line 258 | void SimSetup::createSim( void ){
258      }
259      simnfo->box_z = the_globals->getBoxZ();
260    }
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
261  
223 –    delete fileInit;    
224 –  }
225 –  else{
226 –    initFromBass();
227 –  }
262  
263 <  if( the_globals->haveFinalConfig() ){
264 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
265 <  }
232 <  else{
233 <    strcpy( simnfo->finalName, inFileName );
234 <    char* endTest;
235 <    int nameLength = strlen( simnfo->finalName );
236 <    endTest = &(simnfo->finalName[nameLength - 5]);
237 <    if( !strcmp( endTest, ".bass" ) ){
238 <      strcpy( endTest, ".eor" );
239 <    }
240 <    else if( !strcmp( endTest, ".BASS" ) ){
241 <      strcpy( endTest, ".eor" );
242 <    }
243 <    else{
244 <      endTest = &(simnfo->finalName[nameLength - 4]);
245 <      if( !strcmp( endTest, ".bss" ) ){
246 <        strcpy( endTest, ".eor" );
247 <      }
248 <      else if( !strcmp( endTest, ".mdl" ) ){
249 <        strcpy( endTest, ".eor" );
250 <      }
251 <      else{
252 <        strcat( simnfo->finalName, ".eor" );
253 <      }
254 <    }
255 <  }
256 <    
257 <  // make the sample and status out names
263 > //   if( the_globals->haveInitialConfig() ){
264 > //        InitializeFromFile* fileInit;
265 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
266  
267 <  strcpy( simnfo->sampleName, inFileName );
268 <  char* endTest;
269 <  int nameLength = strlen( simnfo->sampleName );
270 <  endTest = &(simnfo->sampleName[nameLength - 5]);
271 <  if( !strcmp( endTest, ".bass" ) ){
272 <    strcpy( endTest, ".dump" );
273 <  }
274 <  else if( !strcmp( endTest, ".BASS" ) ){
275 <    strcpy( endTest, ".dump" );
276 <  }
277 <  else{
278 <    endTest = &(simnfo->sampleName[nameLength - 4]);
279 <    if( !strcmp( endTest, ".bss" ) ){
280 <      strcpy( endTest, ".dump" );
281 <    }
282 <    else if( !strcmp( endTest, ".mdl" ) ){
283 <      strcpy( endTest, ".dump" );
284 <    }
285 <    else{
286 <      strcat( simnfo->sampleName, ".dump" );
287 <    }
288 <  }
289 <  
290 <  strcpy( simnfo->statusName, inFileName );
291 <  nameLength = strlen( simnfo->statusName );
292 <  endTest = &(simnfo->statusName[nameLength - 5]);
293 <  if( !strcmp( endTest, ".bass" ) ){
294 <    strcpy( endTest, ".stat" );
295 <  }
296 <  else if( !strcmp( endTest, ".BASS" ) ){
297 <    strcpy( endTest, ".stat" );
298 <  }
299 <  else{
300 <    endTest = &(simnfo->statusName[nameLength - 4]);
301 <    if( !strcmp( endTest, ".bss" ) ){
302 <      strcpy( endTest, ".stat" );
303 <    }
304 <    else if( !strcmp( endTest, ".mdl" ) ){
305 <      strcpy( endTest, ".stat" );
306 <    }
307 <    else{
308 <      strcat( simnfo->statusName, ".stat" );
309 <    }
310 <  }
311 <  
267 > //     fileInit->read_xyz( simnfo ); // default velocities on
268 >
269 > //     delete fileInit;
270 > //   }
271 > //   else{
272 >
273 >    initFromBass();
274 >
275 >
276 > //   }
277 >
278 > //   if( the_globals->haveFinalConfig() ){
279 > //     strcpy( simnfo->finalName, the_globals->getFinalConfig() );
280 > //   }
281 > //   else{
282 > //     strcpy( simnfo->finalName, inFileName );
283 > //     char* endTest;
284 > //     int nameLength = strlen( simnfo->finalName );
285 > //     endTest = &(simnfo->finalName[nameLength - 5]);
286 > //     if( !strcmp( endTest, ".bass" ) ){
287 > //       strcpy( endTest, ".eor" );
288 > //     }
289 > //     else if( !strcmp( endTest, ".BASS" ) ){
290 > //       strcpy( endTest, ".eor" );
291 > //     }
292 > //     else{
293 > //       endTest = &(simnfo->finalName[nameLength - 4]);
294 > //       if( !strcmp( endTest, ".bss" ) ){
295 > //      strcpy( endTest, ".eor" );
296 > //       }
297 > //       else if( !strcmp( endTest, ".mdl" ) ){
298 > //      strcpy( endTest, ".eor" );
299 > //       }
300 > //       else{
301 > //      strcat( simnfo->finalName, ".eor" );
302 > //       }
303 > //     }
304 > //   }
305 >
306 > //   // make the sample and status out names
307 >
308 > //   strcpy( simnfo->sampleName, inFileName );
309 > //   char* endTest;
310 > //   int nameLength = strlen( simnfo->sampleName );
311 > //   endTest = &(simnfo->sampleName[nameLength - 5]);
312 > //   if( !strcmp( endTest, ".bass" ) ){
313 > //     strcpy( endTest, ".dump" );
314 > //   }
315 > //   else if( !strcmp( endTest, ".BASS" ) ){
316 > //     strcpy( endTest, ".dump" );
317 > //   }
318 > //   else{
319 > //     endTest = &(simnfo->sampleName[nameLength - 4]);
320 > //     if( !strcmp( endTest, ".bss" ) ){
321 > //       strcpy( endTest, ".dump" );
322 > //     }
323 > //     else if( !strcmp( endTest, ".mdl" ) ){
324 > //       strcpy( endTest, ".dump" );
325 > //     }
326 > //     else{
327 > //       strcat( simnfo->sampleName, ".dump" );
328 > //     }
329 > //   }
330 >
331 > //   strcpy( simnfo->statusName, inFileName );
332 > //   nameLength = strlen( simnfo->statusName );
333 > //   endTest = &(simnfo->statusName[nameLength - 5]);
334 > //   if( !strcmp( endTest, ".bass" ) ){
335 > //     strcpy( endTest, ".stat" );
336 > //   }
337 > //   else if( !strcmp( endTest, ".BASS" ) ){
338 > //     strcpy( endTest, ".stat" );
339 > //   }
340 > //   else{
341 > //     endTest = &(simnfo->statusName[nameLength - 4]);
342 > //     if( !strcmp( endTest, ".bss" ) ){
343 > //       strcpy( endTest, ".stat" );
344 > //     }
345 > //     else if( !strcmp( endTest, ".mdl" ) ){
346 > //       strcpy( endTest, ".stat" );
347 > //     }
348 > //     else{
349 > //       strcat( simnfo->statusName, ".stat" );
350 > //     }
351 > //   }
352 >
353 >
354    // set the status, sample, and themal kick times
355  
356    if( the_globals->haveSampleTime() ){
357 <    simnfo->sampleTime = the_globals->getSampleTime();
357 >    simnfo->sampleTime = the_globals->getSampleTime();
358      simnfo->statusTime = simnfo->sampleTime;
359      simnfo->thermalTime = simnfo->sampleTime;
360    }
361    else{
362 <    simnfo->sampleTime = the_globals->getRunTime();
362 >    simnfo->sampleTime = the_globals->getRunTime();
363      simnfo->statusTime = simnfo->sampleTime;
364      simnfo->thermalTime = simnfo->sampleTime;
365    }
# Line 325 | Line 375 | void SimSetup::createSim( void ){
375    // check for the temperature set flag
376  
377    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
378 <  
379 <    
378 >
379 >
380    // make the longe range forces and the integrator
381 <  
381 >
382    new AllLong( simnfo );
383 <      
383 >
384    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
385    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
386    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
387   }
388  
389   void SimSetup::makeAtoms( void ){
390 <  
390 >
391    int i, j, k, index;
392    double ux, uy, uz, uSqr, u;
393    AtomStamp* current_atom;
394    DirectionalAtom* dAtom;
395 +  int molIndex, molStart, molEnd, nMemb;
396  
397 +
398 +  molIndex = 0;
399    index = 0;
400    for( i=0; i<n_components; i++ ){
401 <    
401 >
402      for( j=0; j<components_nmol[i]; j++ ){
403 <      
403 >
404 >      molStart = index;
405 >      nMemb = comp_stamps[i]->getNAtoms();
406        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
407 <        
407 >
408          current_atom = comp_stamps[i]->getAtom( k );
409 <        if( current_atom->haveOrientation() ){
409 >        if( current_atom->haveOrientation() ){
410  
411 <          dAtom = new DirectionalAtom;
411 >          dAtom = new DirectionalAtom(index);
412            simnfo->n_oriented++;
413            the_atoms[index] = dAtom;
414 <      
414 >
415            ux = current_atom->getOrntX();
416            uy = current_atom->getOrntY();
417            uz = current_atom->getOrntZ();
418 <          
418 >
419            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
420 <          
420 >
421            u = sqrt( uSqr );
422            ux = ux / u;
423            uy = uy / u;
424            uz = uz / u;
425 <          
425 >
426            dAtom->setSUx( ux );
427            dAtom->setSUy( uy );
428            dAtom->setSUz( uz );
429          }
430          else{
431 <          the_atoms[index] = new GeneralAtom;
431 >          the_atoms[index] = new GeneralAtom(index);
432          }
433          the_atoms[index]->setType( current_atom->getType() );
434          the_atoms[index]->setIndex( index );
435 <        
435 >
436          // increment the index and repeat;
437          index++;
438        }
439 +
440 +      molEnd = index -1;
441 +      the_molecules[molIndex].setNMembers( nMemb );
442 +      the_molecules[molIndex].setStartAtom( molStart );
443 +      the_molecules[molIndex].setEndAtom( molEnd );
444 +      molIndex++;
445 +
446      }
447    }
448 <  
448 >
449    the_ff->initializeAtoms();
450   }
451  
# Line 397 | Line 459 | void SimSetup::makeBonds( void ){
459    index = 0;
460    offset = 0;
461    for( i=0; i<n_components; i++ ){
462 <    
462 >
463      for( j=0; j<components_nmol[i]; j++ ){
464 <      
464 >
465        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
466 <        
466 >
467          current_bond = comp_stamps[i]->getBond( k );
468          the_bonds[index].a = current_bond->getA() + offset;
469          the_bonds[index].b = current_bond->getB() + offset;
# Line 415 | Line 477 | void SimSetup::makeBonds( void ){
477        offset += comp_stamps[i]->getNAtoms();
478      }
479    }
480 <  
480 >
481    the_ff->initializeBonds( the_bonds );
482   }
483  
# Line 429 | Line 491 | void SimSetup::makeBends( void ){
491    index = 0;
492    offset = 0;
493    for( i=0; i<n_components; i++ ){
494 <    
494 >
495      for( j=0; j<components_nmol[i]; j++ ){
496 <      
496 >
497        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
498 <        
498 >
499          current_bend = comp_stamps[i]->getBend( k );
500          the_bends[index].a = current_bend->getA() + offset;
501          the_bends[index].b = current_bend->getB() + offset;
# Line 448 | Line 510 | void SimSetup::makeBends( void ){
510        offset += comp_stamps[i]->getNAtoms();
511      }
512    }
513 <  
513 >
514    the_ff->initializeBends( the_bends );
515   }
516  
# Line 462 | Line 524 | void SimSetup::makeTorsions( void ){
524    index = 0;
525    offset = 0;
526    for( i=0; i<n_components; i++ ){
527 <    
527 >
528      for( j=0; j<components_nmol[i]; j++ ){
529 <      
529 >
530        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
531 <        
531 >
532          current_torsion = comp_stamps[i]->getTorsion( k );
533          the_torsions[index].a = current_torsion->getA() + offset;
534          the_torsions[index].b = current_torsion->getB() + offset;
# Line 482 | Line 544 | void SimSetup::makeTorsions( void ){
544        offset += comp_stamps[i]->getNAtoms();
545      }
546    }
547 <  
547 >
548    the_ff->initializeTorsions( the_torsions );
549   }
550  
489 – void SimSetup::makeMolecules( void ){
490 –
491 –  //empy for now
492 – }
493 –
551   void SimSetup::initFromBass( void ){
552  
553    int i, j, k;
# Line 528 | Line 585 | void SimSetup::initFromBass( void ){
585      celly = simnfo->box_y / temp3;
586      cellz = simnfo->box_z / temp3;
587    }
588 <  
588 >
589    current_mol = 0;
590    current_comp_mol = 0;
591    current_comp = 0;
592    current_atom_ndx = 0;
593 <  
593 >
594    for( i=0; i < n_cells ; i++ ){
595      for( j=0; j < n_cells; j++ ){
596        for( k=0; k < n_cells; k++ ){
597 <        
597 >
598          makeElement( i * cellx,
599                       j * celly,
600                       k * cellz );
601 <        
601 >
602          makeElement( i * cellx + 0.5 * cellx,
603                       j * celly + 0.5 * celly,
604                       k * cellz );
605 <        
605 >
606          makeElement( i * cellx,
607                       j * celly + 0.5 * celly,
608                       k * cellz + 0.5 * cellz );
609 <        
609 >
610          makeElement( i * cellx + 0.5 * cellx,
611                       j * celly,
612                       k * cellz + 0.5 * cellz );
# Line 559 | Line 616 | void SimSetup::initFromBass( void ){
616  
617    if( have_extra ){
618      done = 0;
619 <    
619 >
620      int start_ndx;
621      for( i=0; i < (n_cells+1) && !done; i++ ){
622        for( j=0; j < (n_cells+1) && !done; j++ ){
623 <        
623 >
624          if( i < n_cells ){
625 <          
625 >
626            if( j < n_cells ){
627              start_ndx = n_cells;
628            }
629            else start_ndx = 0;
630          }
631          else start_ndx = 0;
632 <        
632 >
633          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
634 <          
634 >
635            makeElement( i * cellx,
636                         j * celly,
637                         k * cellz );
638            done = ( current_mol >= tot_nmol );
639 <          
639 >
640            if( !done && n_per_extra > 1 ){
641              makeElement( i * cellx + 0.5 * cellx,
642                           j * celly + 0.5 * celly,
643                           k * cellz );
644              done = ( current_mol >= tot_nmol );
645            }
646 <          
646 >
647            if( !done && n_per_extra > 2){
648              makeElement( i * cellx,
649                           j * celly + 0.5 * celly,
650                           k * cellz + 0.5 * cellz );
651              done = ( current_mol >= tot_nmol );
652            }
653 <          
653 >
654            if( !done && n_per_extra > 3){
655              makeElement( i * cellx + 0.5 * cellx,
656                           j * celly,
# Line 604 | Line 661 | void SimSetup::initFromBass( void ){
661        }
662      }
663    }
664 <  
665 <  
664 >
665 >
666    for( i=0; i<simnfo->n_atoms; i++ ){
667      simnfo->atoms[i]->set_vx( 0.0 );
668      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 678 | void SimSetup::makeElement( double x, double y, double
678    double rotMat[3][3];
679  
680    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
681 <    
681 >
682      current_atom = comp_stamps[current_comp]->getAtom( k );
683      if( !current_atom->havePosition() ){
684        std::cerr << "Component " << comp_stamps[current_comp]->getID()
# Line 631 | Line 688 | void SimSetup::makeElement( double x, double y, double
688                  << " position.\n";
689        exit(8);
690      }
691 <    
691 >
692      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
693      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
694      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
695 <    
695 >
696      if( the_atoms[current_atom_ndx]->isDirectional() ){
697 <      
697 >
698        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
699 <      
699 >
700        rotMat[0][0] = 1.0;
701        rotMat[0][1] = 0.0;
702        rotMat[0][2] = 0.0;
# Line 657 | Line 714 | void SimSetup::makeElement( double x, double y, double
714  
715      current_atom_ndx++;
716    }
717 <  
717 >
718    current_mol++;
719    current_comp_mol++;
720  
721    if( current_comp_mol >= components_nmol[current_comp] ){
722 <    
722 >
723      current_comp_mol = 0;
724      current_comp++;
725    }

Diff Legend

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