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 128 by chuckv, Thu Oct 3 21:52:46 2002 UTC

# Line 6 | Line 6
6   #include "parse_me.h"
7   #include "LRI.hpp"
8   #include "Integrator.hpp"
9 + #include "mpiInterface.h"
10  
11   SimSetup::SimSetup(){
12    stamps = new MakeStamps();
# Line 21 | Line 22 | void SimSetup::parseFile( char* fileName ){
22  
23    inFileName = fileName;
24    set_interface_stamps( stamps, globals );
25 + #ifdef MPI
26 +  mpiEventInit();
27 + #endif
28    yacc_BASS( fileName );
29 + #ifdef MPI
30 +  throwMPIEvent(NULL);
31 + #endif
32 +
33   }
34  
35 + #ifdef MPI
36 + void SimSetup::receiveParse(void){
37 +
38 +    set_interface_stamps( stamps, globals );
39 +    mpiEventInit();
40 +    mpiEventLoop();
41 +
42 + }
43 + #endif
44 +
45 + void SimSetup::testMe(void){
46 +  bassDiag* dumpMe = new bassDiag(globals,stamps);
47 +  dumpMe->dumpStamps();
48 +  delete dumpMe;
49 + }
50 +
51   void SimSetup::createSim( void ){
52  
53    MakeStamps *the_stamps;
# Line 43 | Line 67 | void SimSetup::createSim( void ){
67    n_components = the_globals->getNComponents();
68    strcpy( force_field, the_globals->getForceField() );
69    strcpy( ensemble, the_globals->getEnsemble() );
70 <  
70 >
71    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
72    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
73    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
74    else{
75 <    std::cerr<< "SimSetup Error. Unrecognized force field -> "
75 >    std::cerr<< "SimSetup Error. Unrecognized force field -> "
76               << force_field << "\n";
77      exit(8);
78    }
# Line 57 | Line 81 | void SimSetup::createSim( void ){
81    the_components = the_globals->getComponents();
82    components_nmol = new int[n_components];
83    comp_stamps = new MoleculeStamp*[n_components];
84 <  
84 >
85    if( !the_globals->haveNMol() ){
86 <    // we don't have the total number of molecules, so we assume it is
86 >    // we don't have the total number of molecules, so we assume it is
87      // given in each component
88  
89      tot_nmol = 0;
90      for( i=0; i<n_components; i++ ){
91 <      
91 >
92        if( !the_components[i]->haveNMol() ){
93          // we have a problem
94          std::cerr << "SimSetup Error. No global NMol or component NMol"
# Line 78 | Line 102 | void SimSetup::createSim( void ){
102    }
103    else{
104      std::cerr << "NOT A SUPPORTED FEATURE\n";
105 <    
105 >
106   //     tot_nmol = the_globals->getNMol();
107 <    
107 >
108   //     //we have the total number of molecules, now we check for molfractions
109   //     for( i=0; i<n_components; i++ ){
110 <      
110 >
111   //       if( !the_components[i]->haveMolFraction() ){
112 <        
112 >
113   //      if( !the_components[i]->haveNMol() ){
114   //        //we have a problem
115   //        std::cerr << "SimSetup error. Neither molFraction nor "
# Line 97 | Line 121 | void SimSetup::createSim( void ){
121  
122    for( i=0; i<n_components; i++ ){
123  
124 <    comp_stamps[i] =
124 >    comp_stamps[i] =
125        the_stamps->getMolecule( the_components[i]->getType() );
126    }
127  
104 –  
128  
129 +
130    // caclulate the number of atoms, bonds, bends and torsions
131  
132    tot_atoms = 0;
# Line 110 | Line 134 | void SimSetup::createSim( void ){
134    tot_bends = 0;
135    tot_torsions = 0;
136    for( i=0; i<n_components; i++ ){
137 <    
137 >
138      tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
139      tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
140      tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
141      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
142    }
143 <  
143 >
144    tot_SRI = tot_bonds + tot_bends + tot_torsions;
145 <  
145 >
146    simnfo->n_atoms = tot_atoms;
147    simnfo->n_bonds = tot_bonds;
148    simnfo->n_bends = tot_bends;
# Line 126 | Line 150 | void SimSetup::createSim( void ){
150    simnfo->n_SRI = tot_SRI;
151  
152    // create the atom and short range interaction arrays
153 <  
153 >
154    the_atoms = new Atom*[tot_atoms];
155 <  //  the_molecules = new Molecule[tot_nmol];
156 <  
157 <  
155 >  the_molecules = new Molecule[tot_nmol];
156 >
157 >
158    if( tot_SRI ){
159      the_sris = new SRI*[tot_SRI];
160      the_excludes = new ex_pair[tot_SRI];
# Line 143 | Line 167 | void SimSetup::createSim( void ){
167    simnfo->n_exclude = tot_SRI;
168    simnfo->excludes = the_excludes;
169  
170 +
171    // initialize the arrays
172 <  
172 >
173    the_ff->setSimInfo( simnfo );
174 <    
174 >
175    makeAtoms();
176  
177    if( tot_bonds ){
# Line 187 | Line 212 | void SimSetup::createSim( void ){
212      simnfo->box_z = the_globals->getBox();
213    }
214    else if( the_globals->haveDensity() ){
215 <    
215 >
216      double vol;
217      vol = (double)tot_nmol / the_globals->getDensity();
218      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 213 | Line 238 | void SimSetup::createSim( void ){
238      }
239      simnfo->box_z = the_globals->getBoxZ();
240    }
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
241  
242 <    delete fileInit;    
243 <  }
244 <  else{
242 >
243 > //   if( the_globals->haveInitialConfig() ){
244 > //        InitializeFromFile* fileInit;
245 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
246 >
247 > //     fileInit->read_xyz( simnfo ); // default velocities on
248 >
249 > //     delete fileInit;
250 > //   }
251 > //   else{
252 >
253      initFromBass();
227 –  }
254  
255 <  if( the_globals->haveFinalConfig() ){
256 <    strcpy( simnfo->finalName, the_globals->getFinalConfig() );
231 <  }
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
255 >
256 > //   }
257  
258 <  strcpy( simnfo->sampleName, inFileName );
259 <  char* endTest;
260 <  int nameLength = strlen( simnfo->sampleName );
261 <  endTest = &(simnfo->sampleName[nameLength - 5]);
262 <  if( !strcmp( endTest, ".bass" ) ){
263 <    strcpy( endTest, ".dump" );
264 <  }
265 <  else if( !strcmp( endTest, ".BASS" ) ){
266 <    strcpy( endTest, ".dump" );
267 <  }
268 <  else{
269 <    endTest = &(simnfo->sampleName[nameLength - 4]);
270 <    if( !strcmp( endTest, ".bss" ) ){
271 <      strcpy( endTest, ".dump" );
272 <    }
273 <    else if( !strcmp( endTest, ".mdl" ) ){
274 <      strcpy( endTest, ".dump" );
275 <    }
276 <    else{
277 <      strcat( simnfo->sampleName, ".dump" );
278 <    }
279 <  }
280 <  
281 <  strcpy( simnfo->statusName, inFileName );
282 <  nameLength = strlen( simnfo->statusName );
283 <  endTest = &(simnfo->statusName[nameLength - 5]);
284 <  if( !strcmp( endTest, ".bass" ) ){
285 <    strcpy( endTest, ".stat" );
286 <  }
287 <  else if( !strcmp( endTest, ".BASS" ) ){
288 <    strcpy( endTest, ".stat" );
289 <  }
290 <  else{
291 <    endTest = &(simnfo->statusName[nameLength - 4]);
292 <    if( !strcmp( endTest, ".bss" ) ){
293 <      strcpy( endTest, ".stat" );
294 <    }
295 <    else if( !strcmp( endTest, ".mdl" ) ){
296 <      strcpy( endTest, ".stat" );
297 <    }
298 <    else{
299 <      strcat( simnfo->statusName, ".stat" );
300 <    }
301 <  }
302 <  
258 > //   if( the_globals->haveFinalConfig() ){
259 > //     strcpy( simnfo->finalName, the_globals->getFinalConfig() );
260 > //   }
261 > //   else{
262 > //     strcpy( simnfo->finalName, inFileName );
263 > //     char* endTest;
264 > //     int nameLength = strlen( simnfo->finalName );
265 > //     endTest = &(simnfo->finalName[nameLength - 5]);
266 > //     if( !strcmp( endTest, ".bass" ) ){
267 > //       strcpy( endTest, ".eor" );
268 > //     }
269 > //     else if( !strcmp( endTest, ".BASS" ) ){
270 > //       strcpy( endTest, ".eor" );
271 > //     }
272 > //     else{
273 > //       endTest = &(simnfo->finalName[nameLength - 4]);
274 > //       if( !strcmp( endTest, ".bss" ) ){
275 > //      strcpy( endTest, ".eor" );
276 > //       }
277 > //       else if( !strcmp( endTest, ".mdl" ) ){
278 > //      strcpy( endTest, ".eor" );
279 > //       }
280 > //       else{
281 > //      strcat( simnfo->finalName, ".eor" );
282 > //       }
283 > //     }
284 > //   }
285 >
286 > //   // make the sample and status out names
287 >
288 > //   strcpy( simnfo->sampleName, inFileName );
289 > //   char* endTest;
290 > //   int nameLength = strlen( simnfo->sampleName );
291 > //   endTest = &(simnfo->sampleName[nameLength - 5]);
292 > //   if( !strcmp( endTest, ".bass" ) ){
293 > //     strcpy( endTest, ".dump" );
294 > //   }
295 > //   else if( !strcmp( endTest, ".BASS" ) ){
296 > //     strcpy( endTest, ".dump" );
297 > //   }
298 > //   else{
299 > //     endTest = &(simnfo->sampleName[nameLength - 4]);
300 > //     if( !strcmp( endTest, ".bss" ) ){
301 > //       strcpy( endTest, ".dump" );
302 > //     }
303 > //     else if( !strcmp( endTest, ".mdl" ) ){
304 > //       strcpy( endTest, ".dump" );
305 > //     }
306 > //     else{
307 > //       strcat( simnfo->sampleName, ".dump" );
308 > //     }
309 > //   }
310 >
311 > //   strcpy( simnfo->statusName, inFileName );
312 > //   nameLength = strlen( simnfo->statusName );
313 > //   endTest = &(simnfo->statusName[nameLength - 5]);
314 > //   if( !strcmp( endTest, ".bass" ) ){
315 > //     strcpy( endTest, ".stat" );
316 > //   }
317 > //   else if( !strcmp( endTest, ".BASS" ) ){
318 > //     strcpy( endTest, ".stat" );
319 > //   }
320 > //   else{
321 > //     endTest = &(simnfo->statusName[nameLength - 4]);
322 > //     if( !strcmp( endTest, ".bss" ) ){
323 > //       strcpy( endTest, ".stat" );
324 > //     }
325 > //     else if( !strcmp( endTest, ".mdl" ) ){
326 > //       strcpy( endTest, ".stat" );
327 > //     }
328 > //     else{
329 > //       strcat( simnfo->statusName, ".stat" );
330 > //     }
331 > //   }
332 >
333 >
334    // set the status, sample, and themal kick times
335  
336    if( the_globals->haveSampleTime() ){
337 <    simnfo->sampleTime = the_globals->getSampleTime();
337 >    simnfo->sampleTime = the_globals->getSampleTime();
338      simnfo->statusTime = simnfo->sampleTime;
339      simnfo->thermalTime = simnfo->sampleTime;
340    }
341    else{
342 <    simnfo->sampleTime = the_globals->getRunTime();
342 >    simnfo->sampleTime = the_globals->getRunTime();
343      simnfo->statusTime = simnfo->sampleTime;
344      simnfo->thermalTime = simnfo->sampleTime;
345    }
# Line 325 | Line 355 | void SimSetup::createSim( void ){
355    // check for the temperature set flag
356  
357    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
358 <  
359 <    
358 >
359 >
360    // make the longe range forces and the integrator
361 <  
361 >
362    new AllLong( simnfo );
363 <      
363 >
364    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
365    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
366    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
367   }
368  
369   void SimSetup::makeAtoms( void ){
370 <  
370 >
371    int i, j, k, index;
372    double ux, uy, uz, uSqr, u;
373    AtomStamp* current_atom;
374    DirectionalAtom* dAtom;
375 +  int molIndex, molStart, molEnd, nMemb;
376  
377 +
378 +  molIndex = 0;
379    index = 0;
380    for( i=0; i<n_components; i++ ){
381 <    
381 >
382      for( j=0; j<components_nmol[i]; j++ ){
383 <      
383 >
384 >      molStart = index;
385 >      nMemb = comp_stamps[i]->getNAtoms();
386        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
387 <        
387 >
388          current_atom = comp_stamps[i]->getAtom( k );
389 <        if( current_atom->haveOrientation() ){
389 >        if( current_atom->haveOrientation() ){
390  
391            dAtom = new DirectionalAtom;
392            simnfo->n_oriented++;
393            the_atoms[index] = dAtom;
394 <      
394 >
395            ux = current_atom->getOrntX();
396            uy = current_atom->getOrntY();
397            uz = current_atom->getOrntZ();
398 <          
398 >
399            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
400 <          
400 >
401            u = sqrt( uSqr );
402            ux = ux / u;
403            uy = uy / u;
404            uz = uz / u;
405 <          
405 >
406            dAtom->setSUx( ux );
407            dAtom->setSUy( uy );
408            dAtom->setSUz( uz );
# Line 377 | Line 412 | void SimSetup::makeAtoms( void ){
412          }
413          the_atoms[index]->setType( current_atom->getType() );
414          the_atoms[index]->setIndex( index );
415 <        
415 >
416          // increment the index and repeat;
417          index++;
418        }
419 +
420 +      molEnd = index -1;
421 +      the_molecules[molIndex].setNMembers( nMemb );
422 +      the_molecules[molIndex].setStartAtom( molStart );
423 +      the_molecules[molIndex].setEndAtom( molEnd );
424 +      molIndex++;
425 +
426      }
427    }
428 <  
428 >
429    the_ff->initializeAtoms();
430   }
431  
# Line 397 | Line 439 | void SimSetup::makeBonds( void ){
439    index = 0;
440    offset = 0;
441    for( i=0; i<n_components; i++ ){
442 <    
442 >
443      for( j=0; j<components_nmol[i]; j++ ){
444 <      
444 >
445        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
446 <        
446 >
447          current_bond = comp_stamps[i]->getBond( k );
448          the_bonds[index].a = current_bond->getA() + offset;
449          the_bonds[index].b = current_bond->getB() + offset;
# Line 415 | Line 457 | void SimSetup::makeBonds( void ){
457        offset += comp_stamps[i]->getNAtoms();
458      }
459    }
460 <  
460 >
461    the_ff->initializeBonds( the_bonds );
462   }
463  
# Line 429 | Line 471 | void SimSetup::makeBends( void ){
471    index = 0;
472    offset = 0;
473    for( i=0; i<n_components; i++ ){
474 <    
474 >
475      for( j=0; j<components_nmol[i]; j++ ){
476 <      
476 >
477        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
478 <        
478 >
479          current_bend = comp_stamps[i]->getBend( k );
480          the_bends[index].a = current_bend->getA() + offset;
481          the_bends[index].b = current_bend->getB() + offset;
# Line 448 | Line 490 | void SimSetup::makeBends( void ){
490        offset += comp_stamps[i]->getNAtoms();
491      }
492    }
493 <  
493 >
494    the_ff->initializeBends( the_bends );
495   }
496  
# Line 462 | Line 504 | void SimSetup::makeTorsions( void ){
504    index = 0;
505    offset = 0;
506    for( i=0; i<n_components; i++ ){
507 <    
507 >
508      for( j=0; j<components_nmol[i]; j++ ){
509 <      
509 >
510        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
511 <        
511 >
512          current_torsion = comp_stamps[i]->getTorsion( k );
513          the_torsions[index].a = current_torsion->getA() + offset;
514          the_torsions[index].b = current_torsion->getB() + offset;
# Line 482 | Line 524 | void SimSetup::makeTorsions( void ){
524        offset += comp_stamps[i]->getNAtoms();
525      }
526    }
527 <  
527 >
528    the_ff->initializeTorsions( the_torsions );
529   }
530  
489 – void SimSetup::makeMolecules( void ){
490 –
491 –  //empy for now
492 – }
493 –
531   void SimSetup::initFromBass( void ){
532  
533    int i, j, k;
# Line 528 | Line 565 | void SimSetup::initFromBass( void ){
565      celly = simnfo->box_y / temp3;
566      cellz = simnfo->box_z / temp3;
567    }
568 <  
568 >
569    current_mol = 0;
570    current_comp_mol = 0;
571    current_comp = 0;
572    current_atom_ndx = 0;
573 <  
573 >
574    for( i=0; i < n_cells ; i++ ){
575      for( j=0; j < n_cells; j++ ){
576        for( k=0; k < n_cells; k++ ){
577 <        
577 >
578          makeElement( i * cellx,
579                       j * celly,
580                       k * cellz );
581 <        
581 >
582          makeElement( i * cellx + 0.5 * cellx,
583                       j * celly + 0.5 * celly,
584                       k * cellz );
585 <        
585 >
586          makeElement( i * cellx,
587                       j * celly + 0.5 * celly,
588                       k * cellz + 0.5 * cellz );
589 <        
589 >
590          makeElement( i * cellx + 0.5 * cellx,
591                       j * celly,
592                       k * cellz + 0.5 * cellz );
# Line 559 | Line 596 | void SimSetup::initFromBass( void ){
596  
597    if( have_extra ){
598      done = 0;
599 <    
599 >
600      int start_ndx;
601      for( i=0; i < (n_cells+1) && !done; i++ ){
602        for( j=0; j < (n_cells+1) && !done; j++ ){
603 <        
603 >
604          if( i < n_cells ){
605 <          
605 >
606            if( j < n_cells ){
607              start_ndx = n_cells;
608            }
609            else start_ndx = 0;
610          }
611          else start_ndx = 0;
612 <        
612 >
613          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
614 <          
614 >
615            makeElement( i * cellx,
616                         j * celly,
617                         k * cellz );
618            done = ( current_mol >= tot_nmol );
619 <          
619 >
620            if( !done && n_per_extra > 1 ){
621              makeElement( i * cellx + 0.5 * cellx,
622                           j * celly + 0.5 * celly,
623                           k * cellz );
624              done = ( current_mol >= tot_nmol );
625            }
626 <          
626 >
627            if( !done && n_per_extra > 2){
628              makeElement( i * cellx,
629                           j * celly + 0.5 * celly,
630                           k * cellz + 0.5 * cellz );
631              done = ( current_mol >= tot_nmol );
632            }
633 <          
633 >
634            if( !done && n_per_extra > 3){
635              makeElement( i * cellx + 0.5 * cellx,
636                           j * celly,
# Line 604 | Line 641 | void SimSetup::initFromBass( void ){
641        }
642      }
643    }
644 <  
645 <  
644 >
645 >
646    for( i=0; i<simnfo->n_atoms; i++ ){
647      simnfo->atoms[i]->set_vx( 0.0 );
648      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 658 | void SimSetup::makeElement( double x, double y, double
658    double rotMat[3][3];
659  
660    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
661 <    
661 >
662      current_atom = comp_stamps[current_comp]->getAtom( k );
663      if( !current_atom->havePosition() ){
664        std::cerr << "Component " << comp_stamps[current_comp]->getID()
# Line 631 | Line 668 | void SimSetup::makeElement( double x, double y, double
668                  << " position.\n";
669        exit(8);
670      }
671 <    
671 >
672      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
673      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
674      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
675 <    
675 >
676      if( the_atoms[current_atom_ndx]->isDirectional() ){
677 <      
677 >
678        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
679 <      
679 >
680        rotMat[0][0] = 1.0;
681        rotMat[0][1] = 0.0;
682        rotMat[0][2] = 0.0;
# Line 657 | Line 694 | void SimSetup::makeElement( double x, double y, double
694  
695      current_atom_ndx++;
696    }
697 <  
697 >
698    current_mol++;
699    current_comp_mol++;
700  
701    if( current_comp_mol >= components_nmol[current_comp] ){
702 <    
702 >
703      current_comp_mol = 0;
704      current_comp++;
705    }

Diff Legend

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