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 145 by mmeineke, Fri Oct 18 15:34:21 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 >  Atom::createArrays(tot_atoms);
160 >  the_molecules = new Molecule[tot_nmol];
161 >
162 >
163    if( tot_SRI ){
164      the_sris = new SRI*[tot_SRI];
165      the_excludes = new ex_pair[tot_SRI];
# Line 143 | Line 172 | void SimSetup::createSim( void ){
172    simnfo->n_exclude = tot_SRI;
173    simnfo->excludes = the_excludes;
174  
175 +
176    // initialize the arrays
177 <  
177 >
178    the_ff->setSimInfo( simnfo );
179 <    
179 >
180    makeAtoms();
181  
182    if( tot_bonds ){
# Line 187 | Line 217 | void SimSetup::createSim( void ){
217      simnfo->box_z = the_globals->getBox();
218    }
219    else if( the_globals->haveDensity() ){
220 <    
220 >
221      double vol;
222      vol = (double)tot_nmol / the_globals->getDensity();
223      simnfo->box_x = pow( vol, ( 1.0 / 3.0 ) );
# Line 213 | Line 243 | void SimSetup::createSim( void ){
243      }
244      simnfo->box_z = the_globals->getBoxZ();
245    }
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
246  
247 <    delete fileInit;    
248 <  }
249 <  else{
247 >
248 > //   if( the_globals->haveInitialConfig() ){
249 > //        InitializeFromFile* fileInit;
250 > //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
251 >
252 > //     fileInit->read_xyz( simnfo ); // default velocities on
253 >
254 > //     delete fileInit;
255 > //   }
256 > //   else{
257 >
258      initFromBass();
227 –  }
259  
229 –  if( the_globals->haveFinalConfig() ){
230 –    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
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 > //   }
262 >
263 > //   if( the_globals->haveFinalConfig() ){
264 > //     strcpy( simnfo->finalName, the_globals->getFinalConfig() );
265 > //   }
266 > //   else{
267 > //     strcpy( simnfo->finalName, inFileName );
268 > //     char* endTest;
269 > //     int nameLength = strlen( simnfo->finalName );
270 > //     endTest = &(simnfo->finalName[nameLength - 5]);
271 > //     if( !strcmp( endTest, ".bass" ) ){
272 > //       strcpy( endTest, ".eor" );
273 > //     }
274 > //     else if( !strcmp( endTest, ".BASS" ) ){
275 > //       strcpy( endTest, ".eor" );
276 > //     }
277 > //     else{
278 > //       endTest = &(simnfo->finalName[nameLength - 4]);
279 > //       if( !strcmp( endTest, ".bss" ) ){
280 > //      strcpy( endTest, ".eor" );
281 > //       }
282 > //       else if( !strcmp( endTest, ".mdl" ) ){
283 > //      strcpy( endTest, ".eor" );
284 > //       }
285 > //       else{
286 > //      strcat( simnfo->finalName, ".eor" );
287 > //       }
288 > //     }
289 > //   }
290 >
291 > //   // make the sample and status out names
292 >
293 > //   strcpy( simnfo->sampleName, inFileName );
294 > //   char* endTest;
295 > //   int nameLength = strlen( simnfo->sampleName );
296 > //   endTest = &(simnfo->sampleName[nameLength - 5]);
297 > //   if( !strcmp( endTest, ".bass" ) ){
298 > //     strcpy( endTest, ".dump" );
299 > //   }
300 > //   else if( !strcmp( endTest, ".BASS" ) ){
301 > //     strcpy( endTest, ".dump" );
302 > //   }
303 > //   else{
304 > //     endTest = &(simnfo->sampleName[nameLength - 4]);
305 > //     if( !strcmp( endTest, ".bss" ) ){
306 > //       strcpy( endTest, ".dump" );
307 > //     }
308 > //     else if( !strcmp( endTest, ".mdl" ) ){
309 > //       strcpy( endTest, ".dump" );
310 > //     }
311 > //     else{
312 > //       strcat( simnfo->sampleName, ".dump" );
313 > //     }
314 > //   }
315 >
316 > //   strcpy( simnfo->statusName, inFileName );
317 > //   nameLength = strlen( simnfo->statusName );
318 > //   endTest = &(simnfo->statusName[nameLength - 5]);
319 > //   if( !strcmp( endTest, ".bass" ) ){
320 > //     strcpy( endTest, ".stat" );
321 > //   }
322 > //   else if( !strcmp( endTest, ".BASS" ) ){
323 > //     strcpy( endTest, ".stat" );
324 > //   }
325 > //   else{
326 > //     endTest = &(simnfo->statusName[nameLength - 4]);
327 > //     if( !strcmp( endTest, ".bss" ) ){
328 > //       strcpy( endTest, ".stat" );
329 > //     }
330 > //     else if( !strcmp( endTest, ".mdl" ) ){
331 > //       strcpy( endTest, ".stat" );
332 > //     }
333 > //     else{
334 > //       strcat( simnfo->statusName, ".stat" );
335 > //     }
336 > //   }
337 >
338 >
339    // set the status, sample, and themal kick times
340  
341    if( the_globals->haveSampleTime() ){
342 <    simnfo->sampleTime = the_globals->getSampleTime();
342 >    simnfo->sampleTime = the_globals->getSampleTime();
343      simnfo->statusTime = simnfo->sampleTime;
344      simnfo->thermalTime = simnfo->sampleTime;
345    }
346    else{
347 <    simnfo->sampleTime = the_globals->getRunTime();
347 >    simnfo->sampleTime = the_globals->getRunTime();
348      simnfo->statusTime = simnfo->sampleTime;
349      simnfo->thermalTime = simnfo->sampleTime;
350    }
# Line 325 | Line 360 | void SimSetup::createSim( void ){
360    // check for the temperature set flag
361  
362    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
363 <  
364 <    
363 >
364 >
365    // make the longe range forces and the integrator
366 <  
366 >
367    new AllLong( simnfo );
368 <      
368 >
369    if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
370    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
371    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
372   }
373  
374   void SimSetup::makeAtoms( void ){
375 <  
375 >
376    int i, j, k, index;
377    double ux, uy, uz, uSqr, u;
378    AtomStamp* current_atom;
379    DirectionalAtom* dAtom;
380 +  int molIndex, molStart, molEnd, nMemb;
381  
382 +
383 +  molIndex = 0;
384    index = 0;
385    for( i=0; i<n_components; i++ ){
386 <    
386 >
387      for( j=0; j<components_nmol[i]; j++ ){
388 <      
388 >
389 >      molStart = index;
390 >      nMemb = comp_stamps[i]->getNAtoms();
391        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
392 <        
392 >
393          current_atom = comp_stamps[i]->getAtom( k );
394 <        if( current_atom->haveOrientation() ){
394 >        if( current_atom->haveOrientation() ){
395  
396 <          dAtom = new DirectionalAtom;
396 >          dAtom = new DirectionalAtom(index);
397            simnfo->n_oriented++;
398            the_atoms[index] = dAtom;
399 <      
399 >
400            ux = current_atom->getOrntX();
401            uy = current_atom->getOrntY();
402            uz = current_atom->getOrntZ();
403 <          
403 >
404            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
405 <          
405 >
406            u = sqrt( uSqr );
407            ux = ux / u;
408            uy = uy / u;
409            uz = uz / u;
410 <          
410 >
411            dAtom->setSUx( ux );
412            dAtom->setSUy( uy );
413            dAtom->setSUz( uz );
414          }
415          else{
416 <          the_atoms[index] = new GeneralAtom;
416 >          the_atoms[index] = new GeneralAtom(index);
417          }
418          the_atoms[index]->setType( current_atom->getType() );
419          the_atoms[index]->setIndex( index );
420 <        
420 >
421          // increment the index and repeat;
422          index++;
423        }
424 +
425 +      molEnd = index -1;
426 +      the_molecules[molIndex].setNMembers( nMemb );
427 +      the_molecules[molIndex].setStartAtom( molStart );
428 +      the_molecules[molIndex].setEndAtom( molEnd );
429 +      molIndex++;
430 +
431      }
432    }
433 <  
433 >
434    the_ff->initializeAtoms();
435   }
436  
# Line 397 | Line 444 | void SimSetup::makeBonds( void ){
444    index = 0;
445    offset = 0;
446    for( i=0; i<n_components; i++ ){
447 <    
447 >
448      for( j=0; j<components_nmol[i]; j++ ){
449 <      
449 >
450        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
451 <        
451 >
452          current_bond = comp_stamps[i]->getBond( k );
453          the_bonds[index].a = current_bond->getA() + offset;
454          the_bonds[index].b = current_bond->getB() + offset;
# Line 415 | Line 462 | void SimSetup::makeBonds( void ){
462        offset += comp_stamps[i]->getNAtoms();
463      }
464    }
465 <  
465 >
466    the_ff->initializeBonds( the_bonds );
467   }
468  
# Line 429 | Line 476 | void SimSetup::makeBends( void ){
476    index = 0;
477    offset = 0;
478    for( i=0; i<n_components; i++ ){
479 <    
479 >
480      for( j=0; j<components_nmol[i]; j++ ){
481 <      
481 >
482        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
483 <        
483 >
484          current_bend = comp_stamps[i]->getBend( k );
485          the_bends[index].a = current_bend->getA() + offset;
486          the_bends[index].b = current_bend->getB() + offset;
# Line 448 | Line 495 | void SimSetup::makeBends( void ){
495        offset += comp_stamps[i]->getNAtoms();
496      }
497    }
498 <  
498 >
499    the_ff->initializeBends( the_bends );
500   }
501  
# Line 462 | Line 509 | void SimSetup::makeTorsions( void ){
509    index = 0;
510    offset = 0;
511    for( i=0; i<n_components; i++ ){
512 <    
512 >
513      for( j=0; j<components_nmol[i]; j++ ){
514 <      
514 >
515        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
516 <        
516 >
517          current_torsion = comp_stamps[i]->getTorsion( k );
518          the_torsions[index].a = current_torsion->getA() + offset;
519          the_torsions[index].b = current_torsion->getB() + offset;
# Line 482 | Line 529 | void SimSetup::makeTorsions( void ){
529        offset += comp_stamps[i]->getNAtoms();
530      }
531    }
532 <  
532 >
533    the_ff->initializeTorsions( the_torsions );
534   }
535  
489 – void SimSetup::makeMolecules( void ){
490 –
491 –  //empy for now
492 – }
493 –
536   void SimSetup::initFromBass( void ){
537  
538    int i, j, k;
# Line 528 | Line 570 | void SimSetup::initFromBass( void ){
570      celly = simnfo->box_y / temp3;
571      cellz = simnfo->box_z / temp3;
572    }
573 <  
573 >
574    current_mol = 0;
575    current_comp_mol = 0;
576    current_comp = 0;
577    current_atom_ndx = 0;
578 <  
578 >
579    for( i=0; i < n_cells ; i++ ){
580      for( j=0; j < n_cells; j++ ){
581        for( k=0; k < n_cells; k++ ){
582 <        
582 >
583          makeElement( i * cellx,
584                       j * celly,
585                       k * cellz );
586 <        
586 >
587          makeElement( i * cellx + 0.5 * cellx,
588                       j * celly + 0.5 * celly,
589                       k * cellz );
590 <        
590 >
591          makeElement( i * cellx,
592                       j * celly + 0.5 * celly,
593                       k * cellz + 0.5 * cellz );
594 <        
594 >
595          makeElement( i * cellx + 0.5 * cellx,
596                       j * celly,
597                       k * cellz + 0.5 * cellz );
# Line 559 | Line 601 | void SimSetup::initFromBass( void ){
601  
602    if( have_extra ){
603      done = 0;
604 <    
604 >
605      int start_ndx;
606      for( i=0; i < (n_cells+1) && !done; i++ ){
607        for( j=0; j < (n_cells+1) && !done; j++ ){
608 <        
608 >
609          if( i < n_cells ){
610 <          
610 >
611            if( j < n_cells ){
612              start_ndx = n_cells;
613            }
614            else start_ndx = 0;
615          }
616          else start_ndx = 0;
617 <        
617 >
618          for( k=start_ndx; k < (n_cells+1) && !done; k++ ){
619 <          
619 >
620            makeElement( i * cellx,
621                         j * celly,
622                         k * cellz );
623            done = ( current_mol >= tot_nmol );
624 <          
624 >
625            if( !done && n_per_extra > 1 ){
626              makeElement( i * cellx + 0.5 * cellx,
627                           j * celly + 0.5 * celly,
628                           k * cellz );
629              done = ( current_mol >= tot_nmol );
630            }
631 <          
631 >
632            if( !done && n_per_extra > 2){
633              makeElement( i * cellx,
634                           j * celly + 0.5 * celly,
635                           k * cellz + 0.5 * cellz );
636              done = ( current_mol >= tot_nmol );
637            }
638 <          
638 >
639            if( !done && n_per_extra > 3){
640              makeElement( i * cellx + 0.5 * cellx,
641                           j * celly,
# Line 604 | Line 646 | void SimSetup::initFromBass( void ){
646        }
647      }
648    }
649 <  
650 <  
649 >
650 >
651    for( i=0; i<simnfo->n_atoms; i++ ){
652      simnfo->atoms[i]->set_vx( 0.0 );
653      simnfo->atoms[i]->set_vy( 0.0 );
# Line 621 | Line 663 | void SimSetup::makeElement( double x, double y, double
663    double rotMat[3][3];
664  
665    for( k=0; k<comp_stamps[current_comp]->getNAtoms(); k++ ){
666 <    
666 >
667      current_atom = comp_stamps[current_comp]->getAtom( k );
668      if( !current_atom->havePosition() ){
669        std::cerr << "Component " << comp_stamps[current_comp]->getID()
# Line 631 | Line 673 | void SimSetup::makeElement( double x, double y, double
673                  << " position.\n";
674        exit(8);
675      }
676 <    
676 >
677      the_atoms[current_atom_ndx]->setX( x + current_atom->getPosX() );
678      the_atoms[current_atom_ndx]->setY( y + current_atom->getPosY() );
679      the_atoms[current_atom_ndx]->setZ( z + current_atom->getPosZ() );
680 <    
680 >
681      if( the_atoms[current_atom_ndx]->isDirectional() ){
682 <      
682 >
683        dAtom = (DirectionalAtom *)the_atoms[current_atom_ndx];
684 <      
684 >
685        rotMat[0][0] = 1.0;
686        rotMat[0][1] = 0.0;
687        rotMat[0][2] = 0.0;
# Line 657 | Line 699 | void SimSetup::makeElement( double x, double y, double
699  
700      current_atom_ndx++;
701    }
702 <  
702 >
703    current_mol++;
704    current_comp_mol++;
705  
706    if( current_comp_mol >= components_nmol[current_comp] ){
707 <    
707 >
708      current_comp_mol = 0;
709      current_comp++;
710    }

Diff Legend

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