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 139 by chuckv, Wed Oct 16 21:07:16 2002 UTC

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

Diff Legend

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