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 195 by chuckv, Thu Dec 5 18:53:40 2002 UTC vs.
Revision 254 by chuckv, Thu Jan 30 20:03:37 2003 UTC

# Line 10 | Line 10
10  
11   #ifdef IS_MPI
12   #include "mpiBASS.h"
13 + #include "mpiSimulation.hpp"
14   #include "bassDiag.hpp"
15   #endif
16  
# Line 73 | Line 74 | void SimSetup::createSim( void ){
74  
75    MakeStamps *the_stamps;
76    Globals* the_globals;
77 <  int i;
77 >  int i, j;
78  
79    // get the stamps and globals;
80    the_stamps = stamps;
# Line 92 | Line 93 | void SimSetup::createSim( void ){
93    if( !strcmp( force_field, "TraPPE" ) ) the_ff = new TraPPEFF();
94    else if( !strcmp( force_field, "DipoleTest" ) ) the_ff = new DipoleTestFF();
95    else if( !strcmp( force_field, "TraPPE_Ex" ) ) the_ff = new TraPPE_ExFF();
96 +  else if( !strcmp( force_field, "LJ" ) ) the_ff = new LJ_FF();
97    else{
98      sprintf( painCave.errMsg,
99               "SimSetup Error. Unrecognized force field -> %s\n",
# Line 161 | Line 163 | void SimSetup::createSim( void ){
163   #endif // is_mpi
164  
165    // make an array of molecule stamps that match the components used.
166 +  // also extract the used stamps out into a separate linked list
167  
168 +  simnfo->nComponents = n_components;
169 +  simnfo->componentsNmol = components_nmol;
170 +  simnfo->compStamps = comp_stamps;
171 +  simnfo->headStamp = new LinkedMolStamp();
172 +  
173 +  char* id;
174 +  LinkedMolStamp* headStamp = simnfo->headStamp;
175 +  LinkedMolStamp* currentStamp = NULL;
176    for( i=0; i<n_components; i++ ){
177  
178 <    comp_stamps[i] =
179 <      the_stamps->getMolecule( the_components[i]->getType() );
178 >    id = the_components[i]->getType();
179 >    comp_stamps[i] = NULL;
180 >    
181 >    // check to make sure the component isn't already in the list
182 >
183 >    comp_stamps[i] = headStamp->match( id );
184 >    if( comp_stamps[i] == NULL ){
185 >      
186 >      // extract the component from the list;
187 >      
188 >      currentStamp = the_stamps->extractMolStamp( id );
189 >      if( currentStamp == NULL ){
190 >        sprintf( painCave.errMsg,
191 >                 "SimSetup error: Component \"%s\" was not found in the "
192 >                 "list of declared molecules\n",
193 >                 id );
194 >        painCave.isFatal = 1;
195 >        simError();
196 >      }
197 >      
198 >      headStamp->add( currentStamp );
199 >      comp_stamps[i] = headStamp->match( id );
200 >    }
201    }
202  
203 + #ifdef IS_MPI
204 +  strcpy( checkPointMsg, "Component stamps successfully extracted\n" );
205 +  MPIcheckPoint();
206 + #endif // is_mpi
207 +  
208  
209  
210 +
211    // caclulate the number of atoms, bonds, bends and torsions
212  
213    tot_atoms = 0;
# Line 177 | Line 215 | void SimSetup::createSim( void ){
215    tot_bends = 0;
216    tot_torsions = 0;
217    for( i=0; i<n_components; i++ ){
218 <
219 <    tot_atoms += components_nmol[i] * comp_stamps[i]->getNAtoms();
220 <    tot_bonds += components_nmol[i] * comp_stamps[i]->getNBonds();
221 <    tot_bends += components_nmol[i] * comp_stamps[i]->getNBends();
218 >    
219 >    tot_atoms +=    components_nmol[i] * comp_stamps[i]->getNAtoms();
220 >    tot_bonds +=    components_nmol[i] * comp_stamps[i]->getNBonds();
221 >    tot_bends +=    components_nmol[i] * comp_stamps[i]->getNBends();
222      tot_torsions += components_nmol[i] * comp_stamps[i]->getNTorsions();
223    }
224  
# Line 198 | Line 236 | void SimSetup::createSim( void ){
236  
237    // divide the molecules among processors here.
238    
239 <  new mpiSimulation( simnfo );
239 >  mpiSim = new mpiSimulation( simnfo );
240    
241 <  simnfo->mpiSim->divideLabor( n_components, comp_stamps, components_nmol );
241 >  fprintf( stderr, "about to call divideLabour.\n" );
242  
243 +  globalIndex = mpiSim->divideLabor();
244 +
245 +  fprintf(stderr, "we're back from divideLabour\n" );
246 +
247 +  // set up the local variables
248 +  
249 +  int localMol, allMol;
250 +  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
251 +  
252 +  allMol = 0;
253 +  localMol = 0;
254 +  local_atoms = 0;
255 +  local_bonds = 0;
256 +  local_bends = 0;
257 +  local_torsions = 0;
258 +  for( i=0; i<n_components; i++ ){
259 +
260 +    for( j=0; j<components_nmol[i]; j++ ){
261 +      
262 +      if( mpiSim->getMyMolStart() <= allMol &&
263 +          allMol <= mpiSim->getMyMolEnd() ){
264 +        
265 +        local_atoms +=    comp_stamps[i]->getNAtoms();
266 +        local_bonds +=    comp_stamps[i]->getNBonds();
267 +        local_bends +=    comp_stamps[i]->getNBends();
268 +        local_torsions += comp_stamps[i]->getNTorsions();
269 +        localMol++;
270 +      }      
271 +      allMol++;
272 +    }
273 +  }
274 +  local_SRI = local_bonds + local_bends + local_torsions;
275 +  
276 +
277 +  simnfo->n_atoms = mpiSim->getMyNlocal();  
278 +  
279 +  if( local_atoms != simnfo->n_atoms ){
280 +    sprintf( painCave.errMsg,
281 +             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
282 +             " localAtom (%d) are note equal.\n",
283 +             simnfo->n_atoms,
284 +             local_atoms );
285 +    painCave.isFatal = 1;
286 +    simError();
287 +  }
288 +
289 +  simnfo->n_bonds = local_bonds;
290 +  simnfo->n_bends = local_bends;
291 +  simnfo->n_torsions = local_torsions;
292 +  simnfo->n_SRI = local_SRI;
293 +  simnfo->n_mol = localMol;
294 +
295 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
296 +  MPIcheckPoint();
297 +  
298 +  
299   #endif // is_mpi
300    
301  
# Line 225 | Line 319 | void SimSetup::createSim( void ){
319    simnfo->excludes = the_excludes;
320  
321  
228 –  // initialize the arrays
229 –
230 –  the_ff->setSimInfo( simnfo );
231 –
232 –  makeAtoms();
233 –
234 –  if( tot_bonds ){
235 –    makeBonds();
236 –  }
237 –
238 –  if( tot_bends ){
239 –    makeBends();
240 –  }
241 –
242 –  if( tot_torsions ){
243 –    makeTorsions();
244 –  }
245 –
246 –
322    // get some of the tricky things that may still be in the globals
323  
324    if( simnfo->n_dipoles ){
# Line 316 | Line 391 | void SimSetup::createSim( void ){
391   #endif // is_mpi
392  
393  
394 +  // initialize the arrays
395  
396 < //   if( the_globals->haveInitialConfig() ){
321 < //        InitializeFromFile* fileInit;
322 < //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
396 >  the_ff->setSimInfo( simnfo );
397  
398 < //     fileInit->read_xyz( simnfo ); // default velocities on
398 >  makeAtoms();
399  
400 < //     delete fileInit;
401 < //   }
402 < //   else{
400 >  if( tot_bonds ){
401 >    makeBonds();
402 >  }
403  
404 +  if( tot_bends ){
405 +    makeBends();
406 +  }
407 +
408 +  if( tot_torsions ){
409 +    makeTorsions();
410 +  }
411 +
412 +
413 +
414 +
415 +
416 +
417 + if( the_globals->haveInitialConfig() ){
418 +
419 +     InitializeFromFile* fileInit;
420 + #ifdef IS_MPI // is_mpi
421 +     if( worldRank == 0 ){
422 + #endif //is_mpi
423 +   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
424 + #ifdef IS_MPI
425 +     }else fileInit = new InitializeFromFile( NULL );
426 + #endif
427 +   fileInit->read_xyz( simnfo ); // default velocities on
428 +
429 +   delete fileInit;
430 + }
431 + else{
432 +
433   #ifdef IS_MPI
434  
435    // no init from bass
# Line 340 | Line 443 | void SimSetup::createSim( void ){
443  
444    initFromBass();
445  
446 < #endif // is_mpi
447 <
446 >
447 > #endif
448 > }
449 >
450   #ifdef IS_MPI
451    strcpy( checkPointMsg, "Successfully read in the initial configuration" );
452    MPIcheckPoint();
# Line 351 | Line 456 | void SimSetup::createSim( void ){
456    
457  
458    
459 <  //   }
459 >
460    
461   #ifdef IS_MPI
462    if( worldRank == 0 ){
# Line 462 | Line 567 | void SimSetup::createSim( void ){
567    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
568  
569  
570 <  // make the longe range forces and the integrator
570 > //   // make the longe range forces and the integrator
571  
572 <  new AllLong( simnfo );
572 > //   new AllLong( simnfo );
573  
574 <  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
574 >  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
575    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
576    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
577 +  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
578 +
579   }
580  
581   void SimSetup::makeAtoms( void ){
# Line 477 | Line 584 | void SimSetup::makeAtoms( void ){
584    double ux, uy, uz, uSqr, u;
585    AtomStamp* current_atom;
586    DirectionalAtom* dAtom;
587 <  int molIndex, molStart, molEnd, nMemb;
587 >  int molIndex, molStart, molEnd, nMemb, lMolIndex;
588  
589 <
589 >  lMolIndex = 0;
590    molIndex = 0;
591    index = 0;
592    for( i=0; i<n_components; i++ ){
593  
594      for( j=0; j<components_nmol[i]; j++ ){
595  
596 <      molStart = index;
597 <      nMemb = comp_stamps[i]->getNAtoms();
598 <      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
596 > #ifdef IS_MPI
597 >      if( mpiSim->getMyMolStart() <= molIndex &&
598 >          molIndex <= mpiSim->getMyMolEnd() ){
599 > #endif // is_mpi        
600  
601 <        current_atom = comp_stamps[i]->getAtom( k );
602 <        if( current_atom->haveOrientation() ){
603 <
604 <          dAtom = new DirectionalAtom(index);
605 <          simnfo->n_oriented++;
606 <          the_atoms[index] = dAtom;
607 <
608 <          ux = current_atom->getOrntX();
609 <          uy = current_atom->getOrntY();
610 <          uz = current_atom->getOrntZ();
611 <
612 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
613 <
614 <          u = sqrt( uSqr );
615 <          ux = ux / u;
616 <          uy = uy / u;
617 <          uz = uz / u;
618 <
619 <          dAtom->setSUx( ux );
620 <          dAtom->setSUy( uy );
621 <          dAtom->setSUz( uz );
601 >        molStart = index;
602 >        nMemb = comp_stamps[i]->getNAtoms();
603 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
604 >          
605 >          current_atom = comp_stamps[i]->getAtom( k );
606 >          if( current_atom->haveOrientation() ){
607 >            
608 >            dAtom = new DirectionalAtom(index);
609 >            simnfo->n_oriented++;
610 >            the_atoms[index] = dAtom;
611 >            
612 >            ux = current_atom->getOrntX();
613 >            uy = current_atom->getOrntY();
614 >            uz = current_atom->getOrntZ();
615 >            
616 >            uSqr = (ux * ux) + (uy * uy) + (uz * uz);
617 >            
618 >            u = sqrt( uSqr );
619 >            ux = ux / u;
620 >            uy = uy / u;
621 >            uz = uz / u;
622 >            
623 >            dAtom->setSUx( ux );
624 >            dAtom->setSUy( uy );
625 >            dAtom->setSUz( uz );
626 >          }
627 >          else{
628 >            the_atoms[index] = new GeneralAtom(index);
629 >          }
630 >          the_atoms[index]->setType( current_atom->getType() );
631 >          the_atoms[index]->setIndex( index );
632 >          
633 >          // increment the index and repeat;
634 >          index++;
635          }
636 <        else{
637 <          the_atoms[index] = new GeneralAtom(index);
638 <        }
639 <        the_atoms[index]->setType( current_atom->getType() );
640 <        the_atoms[index]->setIndex( index );
636 >        
637 >        molEnd = index -1;
638 >        the_molecules[lMolIndex].setNMembers( nMemb );
639 >        the_molecules[lMolIndex].setStartAtom( molStart );
640 >        the_molecules[lMolIndex].setEndAtom( molEnd );
641 >        the_molecules[lMolIndex].setStampID( i );
642 >        lMolIndex++;
643  
644 <        // increment the index and repeat;
522 <        index++;
644 > #ifdef IS_MPI
645        }
646 <
647 <      molEnd = index -1;
526 <      the_molecules[molIndex].setNMembers( nMemb );
527 <      the_molecules[molIndex].setStartAtom( molStart );
528 <      the_molecules[molIndex].setEndAtom( molEnd );
646 > #endif //is_mpi
647 >      
648        molIndex++;
530 –
649      }
650    }
651  
652 + #ifdef IS_MPI
653 +    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
654 +    
655 +    delete[] globalIndex;
656 +
657 +    mpiSim->mpiRefresh();
658 + #endif //IS_MPI
659 +          
660    the_ff->initializeAtoms();
661   }
662  
663   void SimSetup::makeBonds( void ){
664  
665 <  int i, j, k, index, offset;
665 >  int i, j, k, index, offset, molIndex;
666    bond_pair* the_bonds;
667    BondStamp* current_bond;
668  
669    the_bonds = new bond_pair[tot_bonds];
670    index = 0;
671    offset = 0;
672 +  molIndex = 0;
673 +
674    for( i=0; i<n_components; i++ ){
675  
676      for( j=0; j<components_nmol[i]; j++ ){
677  
678 <      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
679 <
680 <        current_bond = comp_stamps[i]->getBond( k );
681 <        the_bonds[index].a = current_bond->getA() + offset;
682 <        the_bonds[index].b = current_bond->getB() + offset;
683 <
684 <        the_excludes[index].i = the_bonds[index].a;
685 <        the_excludes[index].j = the_bonds[index].b;
686 <
687 <        // increment the index and repeat;
688 <        index++;
678 > #ifdef IS_MPI
679 >      if( mpiSim->getMyMolStart() <= molIndex &&
680 >          molIndex <= mpiSim->getMyMolEnd() ){
681 > #endif // is_mpi        
682 >        
683 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
684 >          
685 >          current_bond = comp_stamps[i]->getBond( k );
686 >          the_bonds[index].a = current_bond->getA() + offset;
687 >          the_bonds[index].b = current_bond->getB() + offset;
688 >          
689 >          the_excludes[index].i = the_bonds[index].a;
690 >          the_excludes[index].j = the_bonds[index].b;
691 >          
692 >          // increment the index and repeat;
693 >          index++;
694 >        }
695 >        offset += comp_stamps[i]->getNAtoms();
696 >        
697 > #ifdef IS_MPI
698        }
699 <      offset += comp_stamps[i]->getNAtoms();
700 <    }
699 > #endif //is_mpi
700 >      
701 >      molIndex++;
702 >    }      
703    }
704  
705    the_ff->initializeBonds( the_bonds );
# Line 568 | Line 707 | void SimSetup::makeBends( void ){
707  
708   void SimSetup::makeBends( void ){
709  
710 <  int i, j, k, index, offset;
710 >  int i, j, k, index, offset, molIndex;
711    bend_set* the_bends;
712    BendStamp* current_bend;
713  
714    the_bends = new bend_set[tot_bends];
715    index = 0;
716    offset = 0;
717 +  molIndex = 0;
718    for( i=0; i<n_components; i++ ){
719  
720      for( j=0; j<components_nmol[i]; j++ ){
721  
722 <      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
723 <
724 <        current_bend = comp_stamps[i]->getBend( k );
725 <        the_bends[index].a = current_bend->getA() + offset;
586 <        the_bends[index].b = current_bend->getB() + offset;
587 <        the_bends[index].c = current_bend->getC() + offset;
722 > #ifdef IS_MPI
723 >      if( mpiSim->getMyMolStart() <= molIndex &&
724 >          molIndex <= mpiSim->getMyMolEnd() ){
725 > #endif // is_mpi        
726  
727 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
728 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
729 <
730 <        // increment the index and repeat;
731 <        index++;
727 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
728 >          
729 >          current_bend = comp_stamps[i]->getBend( k );
730 >          the_bends[index].a = current_bend->getA() + offset;
731 >          the_bends[index].b = current_bend->getB() + offset;
732 >          the_bends[index].c = current_bend->getC() + offset;
733 >          
734 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
735 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
736 >          
737 >          // increment the index and repeat;
738 >          index++;
739 >        }
740 >        offset += comp_stamps[i]->getNAtoms();
741 >        
742 > #ifdef IS_MPI
743        }
744 <      offset += comp_stamps[i]->getNAtoms();
744 > #endif //is_mpi
745 >
746 >      molIndex++;
747      }
748    }
749  
# Line 601 | Line 752 | void SimSetup::makeTorsions( void ){
752  
753   void SimSetup::makeTorsions( void ){
754  
755 <  int i, j, k, index, offset;
755 >  int i, j, k, index, offset, molIndex;
756    torsion_set* the_torsions;
757    TorsionStamp* current_torsion;
758  
759    the_torsions = new torsion_set[tot_torsions];
760    index = 0;
761    offset = 0;
762 +  molIndex = 0;
763    for( i=0; i<n_components; i++ ){
764  
765      for( j=0; j<components_nmol[i]; j++ ){
766  
767 + #ifdef IS_MPI
768 +      if( mpiSim->getMyMolStart() <= molIndex &&
769 +          molIndex <= mpiSim->getMyMolEnd() ){
770 + #endif // is_mpi        
771 +
772        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
773  
774          current_torsion = comp_stamps[i]->getTorsion( k );
# Line 627 | Line 784 | void SimSetup::makeTorsions( void ){
784          index++;
785        }
786        offset += comp_stamps[i]->getNAtoms();
787 +
788 + #ifdef IS_MPI
789 +      }
790 + #endif //is_mpi      
791 +
792 +      molIndex++;
793      }
794    }
795  

Diff Legend

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