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 253 by chuckv, Thu Jan 30 15:20:21 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 >  globalIndex = mpiSim->divideLabor();
242 >
243 >  // set up the local variables
244 >  
245 >  int localMol, allMol;
246 >  int local_atoms, local_bonds, local_bends, local_torsions, local_SRI;
247 >  
248 >  allMol = 0;
249 >  localMol = 0;
250 >  local_atoms = 0;
251 >  local_bonds = 0;
252 >  local_bends = 0;
253 >  local_torsions = 0;
254 >  for( i=0; i<n_components; i++ ){
255 >
256 >    for( j=0; j<components_nmol[i]; j++ ){
257 >      
258 >      if( mpiSim->getMyMolStart() <= allMol &&
259 >          allMol <= mpiSim->getMyMolEnd() ){
260 >        
261 >        local_atoms +=    comp_stamps[i]->getNAtoms();
262 >        local_bonds +=    comp_stamps[i]->getNBonds();
263 >        local_bends +=    comp_stamps[i]->getNBends();
264 >        local_torsions += comp_stamps[i]->getNTorsions();
265 >        localMol++;
266 >      }      
267 >      allMol++;
268 >    }
269 >  }
270 >  local_SRI = local_bonds + local_bends + local_torsions;
271 >  
272 >
273 >  simnfo->n_atoms = mpiSim->getMyNlocal();  
274 >  
275 >  if( local_atoms != simnfo->n_atoms ){
276 >    sprintf( painCave.errMsg,
277 >             "SimSetup error: mpiSim's localAtom (%d) and SimSetup's"
278 >             " localAtom (%d) are note equal.\n",
279 >             simnfo->n_atoms,
280 >             local_atoms );
281 >    painCave.isFatal = 1;
282 >    simError();
283 >  }
284 >
285 >  simnfo->n_bonds = local_bonds;
286 >  simnfo->n_bends = local_bends;
287 >  simnfo->n_torsions = local_torsions;
288 >  simnfo->n_SRI = local_SRI;
289 >  simnfo->n_mol = localMol;
290  
291 +  strcpy( checkPointMsg, "Passed nlocal consistency check." );
292 +  MPIcheckPoint();
293 +  
294 +  
295   #endif // is_mpi
296    
297  
# Line 225 | Line 315 | void SimSetup::createSim( void ){
315    simnfo->excludes = the_excludes;
316  
317  
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 –
318    // get some of the tricky things that may still be in the globals
319  
320    if( simnfo->n_dipoles ){
# Line 316 | Line 387 | void SimSetup::createSim( void ){
387   #endif // is_mpi
388  
389  
390 +  // initialize the arrays
391  
392 < //   if( the_globals->haveInitialConfig() ){
321 < //        InitializeFromFile* fileInit;
322 < //     fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
392 >  the_ff->setSimInfo( simnfo );
393  
394 < //     fileInit->read_xyz( simnfo ); // default velocities on
394 >  makeAtoms();
395  
396 < //     delete fileInit;
397 < //   }
398 < //   else{
396 >  if( tot_bonds ){
397 >    makeBonds();
398 >  }
399  
400 +  if( tot_bends ){
401 +    makeBends();
402 +  }
403 +
404 +  if( tot_torsions ){
405 +    makeTorsions();
406 +  }
407 +
408 +
409 +
410 +
411 +
412 +
413 + if( the_globals->haveInitialConfig() ){
414 +
415 +     InitializeFromFile* fileInit;
416 + #ifdef IS_MPI // is_mpi
417 +     if( worldRank == 0 ){
418 + #endif //is_mpi
419 +   fileInit = new InitializeFromFile( the_globals->getInitialConfig() );
420 + #ifdef IS_MPI
421 +     }else fileInit = new InitializeFromFile( NULL );
422 + #endif
423 +   fileInit->read_xyz( simnfo ); // default velocities on
424 +
425 +   delete fileInit;
426 + }
427 + else{
428 +
429   #ifdef IS_MPI
430  
431    // no init from bass
# Line 340 | Line 439 | void SimSetup::createSim( void ){
439  
440    initFromBass();
441  
442 < #endif // is_mpi
443 <
442 >
443 > #endif
444 > }
445 >
446   #ifdef IS_MPI
447    strcpy( checkPointMsg, "Successfully read in the initial configuration" );
448    MPIcheckPoint();
# Line 351 | Line 452 | void SimSetup::createSim( void ){
452    
453  
454    
455 <  //   }
455 >
456    
457   #ifdef IS_MPI
458    if( worldRank == 0 ){
# Line 462 | Line 563 | void SimSetup::createSim( void ){
563    if( the_globals->haveTempSet() ) simnfo->setTemp = the_globals->getTempSet();
564  
565  
566 <  // make the longe range forces and the integrator
566 > //   // make the longe range forces and the integrator
567  
568 <  new AllLong( simnfo );
568 > //   new AllLong( simnfo );
569  
570 <  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo );
570 >  if( !strcmp( force_field, "TraPPE" ) ) new Verlet( *simnfo, the_ff );
571    if( !strcmp( force_field, "DipoleTest" ) ) new Symplectic( simnfo );
572    if( !strcmp( force_field, "TraPPE_Ex" ) ) new Symplectic( simnfo );
573 +  if( !strcmp( force_field, "LJ" ) ) new Verlet( *simnfo, the_ff );
574 +
575   }
576  
577   void SimSetup::makeAtoms( void ){
# Line 477 | Line 580 | void SimSetup::makeAtoms( void ){
580    double ux, uy, uz, uSqr, u;
581    AtomStamp* current_atom;
582    DirectionalAtom* dAtom;
583 <  int molIndex, molStart, molEnd, nMemb;
583 >  int molIndex, molStart, molEnd, nMemb, lMolIndex;
584  
585 <
585 >  lMolIndex = 0;
586    molIndex = 0;
587    index = 0;
588    for( i=0; i<n_components; i++ ){
589  
590      for( j=0; j<components_nmol[i]; j++ ){
591  
592 <      molStart = index;
593 <      nMemb = comp_stamps[i]->getNAtoms();
594 <      for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
592 > #ifdef IS_MPI
593 >      if( mpiSim->getMyMolStart() <= molIndex &&
594 >          molIndex <= mpiSim->getMyMolEnd() ){
595 > #endif // is_mpi        
596  
597 <        current_atom = comp_stamps[i]->getAtom( k );
598 <        if( current_atom->haveOrientation() ){
599 <
600 <          dAtom = new DirectionalAtom(index);
601 <          simnfo->n_oriented++;
602 <          the_atoms[index] = dAtom;
603 <
604 <          ux = current_atom->getOrntX();
605 <          uy = current_atom->getOrntY();
606 <          uz = current_atom->getOrntZ();
607 <
608 <          uSqr = (ux * ux) + (uy * uy) + (uz * uz);
609 <
610 <          u = sqrt( uSqr );
611 <          ux = ux / u;
612 <          uy = uy / u;
613 <          uz = uz / u;
614 <
615 <          dAtom->setSUx( ux );
616 <          dAtom->setSUy( uy );
617 <          dAtom->setSUz( uz );
597 >        molStart = index;
598 >        nMemb = comp_stamps[i]->getNAtoms();
599 >        for( k=0; k<comp_stamps[i]->getNAtoms(); k++ ){
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 );
622 >          }
623 >          else{
624 >            the_atoms[index] = new GeneralAtom(index);
625 >          }
626 >          the_atoms[index]->setType( current_atom->getType() );
627 >          the_atoms[index]->setIndex( index );
628 >          
629 >          // increment the index and repeat;
630 >          index++;
631          }
632 <        else{
633 <          the_atoms[index] = new GeneralAtom(index);
634 <        }
635 <        the_atoms[index]->setType( current_atom->getType() );
636 <        the_atoms[index]->setIndex( index );
632 >        
633 >        molEnd = index -1;
634 >        the_molecules[lMolIndex].setNMembers( nMemb );
635 >        the_molecules[lMolIndex].setStartAtom( molStart );
636 >        the_molecules[lMolIndex].setEndAtom( molEnd );
637 >        the_molecules[lMolIndex].setStampID( i );
638 >        lMolIndex++;
639  
640 <        // increment the index and repeat;
522 <        index++;
640 > #ifdef IS_MPI
641        }
642 <
643 <      molEnd = index -1;
526 <      the_molecules[molIndex].setNMembers( nMemb );
527 <      the_molecules[molIndex].setStartAtom( molStart );
528 <      the_molecules[molIndex].setEndAtom( molEnd );
642 > #endif //is_mpi
643 >      
644        molIndex++;
530 –
645      }
646    }
647  
648 + #ifdef IS_MPI
649 +    for( i=0; i<mpiSim->getMyNlocal(); i++ ) the_atoms[i]->setGlobalIndex( globalIndex[i] );
650 +    
651 +    delete[] globalIndex;
652 + #endif IS_MPI
653 +          
654    the_ff->initializeAtoms();
655   }
656  
657   void SimSetup::makeBonds( void ){
658  
659 <  int i, j, k, index, offset;
659 >  int i, j, k, index, offset, molIndex;
660    bond_pair* the_bonds;
661    BondStamp* current_bond;
662  
663    the_bonds = new bond_pair[tot_bonds];
664    index = 0;
665    offset = 0;
666 +  molIndex = 0;
667 +
668    for( i=0; i<n_components; i++ ){
669  
670      for( j=0; j<components_nmol[i]; j++ ){
549 –
550 –      for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
671  
672 <        current_bond = comp_stamps[i]->getBond( k );
673 <        the_bonds[index].a = current_bond->getA() + offset;
674 <        the_bonds[index].b = current_bond->getB() + offset;
675 <
676 <        the_excludes[index].i = the_bonds[index].a;
677 <        the_excludes[index].j = the_bonds[index].b;
678 <
679 <        // increment the index and repeat;
680 <        index++;
672 > #ifdef IS_MPI
673 >      if( mpiSim->getMyMolStart() <= molIndex &&
674 >          molIndex <= mpiSim->getMyMolEnd() ){
675 > #endif // is_mpi        
676 >        
677 >        for( k=0; k<comp_stamps[i]->getNBonds(); k++ ){
678 >          
679 >          current_bond = comp_stamps[i]->getBond( k );
680 >          the_bonds[index].a = current_bond->getA() + offset;
681 >          the_bonds[index].b = current_bond->getB() + offset;
682 >          
683 >          the_excludes[index].i = the_bonds[index].a;
684 >          the_excludes[index].j = the_bonds[index].b;
685 >          
686 >          // increment the index and repeat;
687 >          index++;
688 >        }
689 >        offset += comp_stamps[i]->getNAtoms();
690 >        
691 > #ifdef IS_MPI
692        }
693 <      offset += comp_stamps[i]->getNAtoms();
694 <    }
693 > #endif is_mpi
694 >      
695 >      molIndex++;
696 >    }      
697    }
698  
699    the_ff->initializeBonds( the_bonds );
# Line 568 | Line 701 | void SimSetup::makeBends( void ){
701  
702   void SimSetup::makeBends( void ){
703  
704 <  int i, j, k, index, offset;
704 >  int i, j, k, index, offset, molIndex;
705    bend_set* the_bends;
706    BendStamp* current_bend;
707  
708    the_bends = new bend_set[tot_bends];
709    index = 0;
710    offset = 0;
711 +  molIndex = 0;
712    for( i=0; i<n_components; i++ ){
713  
714      for( j=0; j<components_nmol[i]; j++ ){
715  
716 <      for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
717 <
718 <        current_bend = comp_stamps[i]->getBend( k );
719 <        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;
716 > #ifdef IS_MPI
717 >      if( mpiSim->getMyMolStart() <= molIndex &&
718 >          molIndex <= mpiSim->getMyMolEnd() ){
719 > #endif // is_mpi        
720  
721 <        the_excludes[index + tot_bonds].i = the_bends[index].a;
722 <        the_excludes[index + tot_bonds].j = the_bends[index].c;
723 <
724 <        // increment the index and repeat;
725 <        index++;
721 >        for( k=0; k<comp_stamps[i]->getNBends(); k++ ){
722 >          
723 >          current_bend = comp_stamps[i]->getBend( k );
724 >          the_bends[index].a = current_bend->getA() + offset;
725 >          the_bends[index].b = current_bend->getB() + offset;
726 >          the_bends[index].c = current_bend->getC() + offset;
727 >          
728 >          the_excludes[index + tot_bonds].i = the_bends[index].a;
729 >          the_excludes[index + tot_bonds].j = the_bends[index].c;
730 >          
731 >          // increment the index and repeat;
732 >          index++;
733 >        }
734 >        offset += comp_stamps[i]->getNAtoms();
735 >        
736 > #ifdef IS_MPI
737        }
738 <      offset += comp_stamps[i]->getNAtoms();
738 > #endif //is_mpi
739 >
740 >      molIndex++;
741      }
742    }
743  
# Line 601 | Line 746 | void SimSetup::makeTorsions( void ){
746  
747   void SimSetup::makeTorsions( void ){
748  
749 <  int i, j, k, index, offset;
749 >  int i, j, k, index, offset, molIndex;
750    torsion_set* the_torsions;
751    TorsionStamp* current_torsion;
752  
753    the_torsions = new torsion_set[tot_torsions];
754    index = 0;
755    offset = 0;
756 +  molIndex = 0;
757    for( i=0; i<n_components; i++ ){
758  
759      for( j=0; j<components_nmol[i]; j++ ){
760  
761 + #ifdef IS_MPI
762 +      if( mpiSim->getMyMolStart() <= molIndex &&
763 +          molIndex <= mpiSim->getMyMolEnd() ){
764 + #endif // is_mpi        
765 +
766        for( k=0; k<comp_stamps[i]->getNTorsions(); k++ ){
767  
768          current_torsion = comp_stamps[i]->getTorsion( k );
# Line 627 | Line 778 | void SimSetup::makeTorsions( void ){
778          index++;
779        }
780        offset += comp_stamps[i]->getNAtoms();
781 +
782 + #ifdef IS_MPI
783 +      }
784 + #endif //is_mpi      
785 +
786 +      molIndex++;
787      }
788    }
789  

Diff Legend

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