ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mdtools/interface_implementation/SimSetup.cpp
(Generate patch)

Comparing:
branches/mmeineke/mdtools/interface_implementation/SimSetup.cpp (file contents), Revision 10 by mmeineke, Tue Jul 9 18:40:59 2002 UTC vs.
trunk/mdtools/interface_implementation/SimSetup.cpp (file contents), Revision 144 by mmeineke, Thu Oct 17 21:59:12 2002 UTC

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

Diff Legend

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