ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-2.0/src/brains/SimInfo.cpp (file contents):
Revision 1722 by tim, Tue Nov 9 23:11:39 2004 UTC vs.
Revision 1735 by tim, Fri Nov 12 17:40:03 2004 UTC

# Line 37 | Line 37 | SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0
37  
38   namespace oopse {
39  
40 < SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
41 <        nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){
40 > SimInfo::SimInfo(const std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
41 >                                ForceField* ff, Globals* globals) :
42 >                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
43 >                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
44 >                                nCutoffGroups_(0), nConstraints_(0), nZConstraint_(0), sman_(NULL) {
45  
46 +    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
47 +    int nCutoffAtoms; // number of atoms belong to cutoff groups
48 +    int ngroups;          //total cutoff groups defined in meta-data file
49 +    MoleculeStamp* molStamp;
50 +    int nMolWithSameStamp;
51 +    CutoffGroupStamp* cgStamp;
52 +    int nAtomsIngroups;
53 +    int nCutoffGroupsInStamp;    
54 +
55 +    nGlobalAtoms_ =  0;
56 +    ngroups = 0;
57 +    
58 +    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
59 +        molStamp = i->first;
60 +        nMolWithSameStamp = i->second;
61 +        
62 +        addMoleculeStamp(molStamp, nMolWithSameStamp);
63 +        
64 +        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
65 +        
66 +        nAtomsIngroups = 0;
67 +        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
68 +        
69 +        for (int j=0; j < nCutoffGroupsInStamp; j++) {
70 +            cgStamp = molStamp->getCutoffGroup(j);
71 +            nAtomsIngroups += cgStamp->getNMembers();
72 +        }
73 +
74 +        ngroups += *nMolWithSameStamp;
75 +        nCutoffAtoms += nAtomsIngroups * nMolWithSameStamp;                
76 +    }
77 +
78 +    //every free atom (atom does not belong to cutoff groups) is a cutoff group
79 +    //therefore the total number of cutoff groups in the system is equal to
80 +    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
81 +    //file plus the number of cutoff groups defined in meta-data file
82 +    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + ngroups;
83 +
84 +    //initialize globalGroupMembership_, every element of this array will be 0
85 +    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
86 +
87 +    nGlobalMols_ = molStampIds_.size();
88 +
89 + #ifdef IS_MPI    
90 +    molToProcMap_.resize(nGlobalMols_);
91 + #endif
92 +    
93   }
94  
95   SimInfo::~SimInfo() {
96 <    MemoryUtils::deleteVectorOfPointer(molecules_);
96 >    //MemoryUtils::deleteVectorOfPointer(molecules_);
97 >
98 >    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
99 >    
100      delete sman_;
101 +    delete globals_;
102 +    delete forceField_;
103  
104   }
105  
106  
107   bool SimInfo::addMolecule(Molecule* mol) {
108 <    std::vector<Molecule*>::iterator i;
109 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
108 >    MoleculeIterator i;
109 >
110 >    i = molecules_.find(mol->getGlobalIndex());
111      if (i != molecules_.end() ) {
56 –        molecules_.push_back(mol);
112  
113 +        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
114 +        
115          nAtoms_ += mol->getNAtoms();
116          nBonds_ += mol->getNBonds();
117          nBends_ += mol->getNBends();
# Line 71 | Line 128 | bool SimInfo::removeMolecule(Molecule* mol) {
128   }
129  
130   bool SimInfo::removeMolecule(Molecule* mol) {
131 <    std::vector<Molecule*>::iterator i;
132 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
131 >    MoleculeIterator i;
132 >    i = molecules_.find(mol->getGlobalIndex());
133  
134      if (i != molecules_.end() ) {
135 <        molecules_.push_back(mol);
135 >
136 >        assert(mol == i->second);
137 >        
138          nAtoms_ -= mol->getNAtoms();
139          nBonds_ -= mol->getNBonds();
140          nBends_ -= mol->getNBends();
# Line 85 | Line 144 | bool SimInfo::removeMolecule(Molecule* mol) {
144          nCutoffGroups_ -= mol->getNCutoffGroups();
145          nConstraints_ -= mol->getNConstraints();
146  
147 +        molecules_.erase(mol->getGlobalIndex());
148 +
149 +        delete mol;
150 +        
151          return true;
152      } else {
153          return false;
# Line 94 | Line 157 | Molecule* SimInfo::beginMolecule(std::vector<Molecule*
157   }    
158  
159          
160 < Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
160 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
161      i = molecules_.begin();
162      return i == molecules_.end() ? NULL : *i;
163   }    
164  
165 < Molecule* SimInfo::nextMolecule(std::vector<Molecule*>::iterator& i) {
165 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
166      ++i;
167      return i == molecules_.end() ? NULL : *i;    
168   }
# Line 107 | Line 170 | void SimInfo::calcNdf() {
170  
171   void SimInfo::calcNdf() {
172      int ndf_local;
173 <    std::vector<Molecule*>::iterator i;
173 >    MoleculeIterator i;
174      std::vector<StuntDouble*>::iterator j;
175      Molecule* mol;
176      StuntDouble* integrableObject;
# Line 140 | Line 203 | void SimInfo::calcNdf() {
203      ndf_ = ndf_local;
204   #endif
205  
206 <    // nZconstraints is global, as are the 3 COM translations for the
206 >    // nZconstraints_ is global, as are the 3 COM translations for the
207      // entire system:
208 <    ndf_ = ndf_ - 3 - nZconstraints;
208 >    ndf_ = ndf_ - 3 - nZconstraints_;
209  
210   }
211  
212   void SimInfo::calcNdfRaw() {
213      int ndfRaw_local;
214  
215 <    std::vector<Molecule*>::iterator i;
215 >    MoleculeIterator i;
216      std::vector<StuntDouble*>::iterator j;
217      Molecule* mol;
218      StuntDouble* integrableObject;
# Line 193 | Line 256 | void SimInfo::calcNdfTrans() {
256      ndfTrans_ = ndfTrans_local;
257   #endif
258  
259 <    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
259 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraints_;
260  
261   }
262  
# Line 282 | Line 345 | void SimInfo::removeExcludePairs(Molecule* mol) {
345          exclude_.removePair(b, c);
346          exclude_.removePair(b, d);
347          exclude_.removePair(c, d);        
348 +    }
349 +
350 + }
351 +
352 +
353 + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
354 +    int curStampId;
355 +
356 +    //index from 0
357 +    curStampId = molStampIds_.size();
358 +
359 +    moleculeStamps_.push_back(molStamp);
360 +    molStampIds_.insert(molStampIds_.end(), nmol, curStampId)
361 + }
362 +
363 + void SimInfo::update() {
364 +
365 + #ifdef IS_MPI
366 +    setupFortranParallel();
367 + #endif
368 +
369 +    setupFortranSim();
370 +
371 +    calcNdf();
372 +    calcNdfRaw();
373 +    calcNdfTrans();
374 + }
375 +
376 + std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
377 +    typename SimInfo::MoleculeIterator mi;
378 +    Molecule* mol;
379 +    typename Molecule::AtomIterator ai;
380 +    Atom* atom;
381 +    std::set<AtomType*> atomTypes;
382 +
383 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
384 +
385 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
386 +            atomTypes.insert(atom->getAtomType());
387 +        }
388 +        
389 +    }
390 +
391 +    return atomTypes;        
392 + }
393 +
394 + void SimInfo::setupSimType() {
395 +    std::set<AtomType*>::iterator i;
396 +    std::set<AtomType*> atomTypes;
397 +    atomTypes = getUniqueAtomTypes();
398 +    
399 +    int useLennardJones = 0;
400 +    int useElectrostatic = 0;
401 +    int useEAM = 0;
402 +    int useCharge = 0;
403 +    int useDirectional = 0;
404 +    int useDipole = 0;
405 +    int useGayBerne = 0;
406 +    int useSticky = 0;
407 +    int useShape = 0;
408 +    int useFLARB = 0; //it is not in AtomType yet
409 +    int useDirectionalAtom = 0;    
410 +    int useElectrostatics = 0;
411 +    //usePBC and useRF are from globals
412 +    bool usePBC = globals_->getPBC();
413 +    bool useRF = globals_->getUseRF();
414 +
415 +    //loop over all of the atom types
416 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
417 +        useLennardJones |= i->isLennardJones();
418 +        useElectrostatic |= i->isElectrostatic();
419 +        useEAM |= i->isEAM();
420 +        useCharge |= i->isCharge();
421 +        useDirectional |= i->isDirectional();
422 +        useDipole |= i->isDipole();
423 +        useGayBerne |= i->isGayBerne();
424 +        useSticky |= i->isSticky();
425 +        useShape |= i->isShape();
426 +    }
427 +
428 +    if (useSticky || useDipole || useGayBerne || useShape) {
429 +        useDirectionalAtom = 1;
430 +    }
431 +
432 +    if (useCharge || useDipole) {
433 +        useElectrostatics = 1;
434 +    }
435 +
436 + #ifdef IS_MPI    
437 +    int temp;
438 +
439 +    temp = usePBC;
440 +    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
441 +
442 +    temp = useDirectionalAtom;
443 +    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
444 +
445 +    temp = useLennardJones;
446 +    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
447 +
448 +    temp = useElectrostatics;
449 +    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
450 +
451 +    temp = useCharge;
452 +    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
453 +
454 +    temp = useDipole;
455 +    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
456 +
457 +    temp = useSticky;
458 +    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
459 +
460 +    temp = useGayBerne;
461 +    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
462 +
463 +    temp = useEAM;
464 +    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
465 +
466 +    temp = useShape;
467 +    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
468 +
469 +    temp = useFLARB;
470 +    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
471 +
472 +    temp = useRF;
473 +    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
474 +    
475 + #endif
476 +
477 +    fInfo_.SIM_uses_PBC = usePBC;    
478 +    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
479 +    fInfo_.SIM_uses_LennardJones = useLennardJones;
480 +    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
481 +    fInfo_.SIM_uses_Charges = useCharge;
482 +    fInfo_.SIM_uses_Dipoles = useDipole;
483 +    fInfo_.SIM_uses_Sticky = useSticky;
484 +    fInfo_.SIM_uses_GayBerne = useGayBerne;
485 +    fInfo_.SIM_uses_EAM = useEAM;
486 +    fInfo_.SIM_uses_Shapes = useShape;
487 +    fInfo_.SIM_uses_FLARB = useFLARB;
488 +    fInfo_.SIM_uses_RF = useRF;
489 +
490 +    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
491 +        fInfo_.dielect = dielectric;
492 +    } else {
493 +        fInfo_.dielect = 0.0;
494      }
495  
496   }
497  
498 + void SimInfo::setupFortranSim() {
499 +    int isError;
500 +    int nExclude;
501 +    std::vector<int> fortranGlobalGroupMembership;
502 +    
503 +    nExclude = exclude_.getSize();
504 +    isError = 0;
505  
506 +    //globalGroupMembership_ is filled by SimCreator    
507 +    for (int i = 0; i < nGlobalAtoms_; i++) {
508 +        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
509 +    }
510 +
511 +    //calculate mass ratio of cutoff group
512 +    std::vector<double> mfact;
513 +    typename SimInfo::MoleculeIterator mi;
514 +    Molecule* mol;
515 +    typename Molecule::CutoffGroupIterator ci;
516 +    CutoffGroup* cg;
517 +    typename Molecule::AtomIterator ai;
518 +    Atom* atom;
519 +    double totalMass;
520 +
521 +    //to avoid memory reallocation, reserve enough space for mfact
522 +    mfact.reserve(getNCutoffGroups());
523 +    
524 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
525 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
526 +
527 +            totalMass = cg->getMass();
528 +            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
529 +                        mfact.push_back(atom->getMass()/totalMass);
530 +            }
531 +
532 +        }      
533 +    }
534 +
535 +    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
536 +    std::vector<int> identArray;
537 +
538 +    //to avoid memory reallocation, reserve enough space identArray
539 +    identArray.reserve(getNAtoms());
540 +    
541 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
542 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
543 +            identArray.push_back(atom->getIdent());
544 +        }
545 +    }    
546 +
547 +    //fill molMembershipArray
548 +    //molMembershipArray is filled by SimCreator    
549 +    std::vector<int> molMembershipArray(nGlobalAtoms_);
550 +    for (int i = 0; i < nGlobalAtoms_; i++) {
551 +        molMembershipArray.push_back(globalMolMembership_[i] + 1);
552 +    }
553 +    
554 +    //setup fortran simulation
555 +    //gloalExcludes and molMembershipArray should go away (They are never used)
556 +    //why the hell fortran need to know molecule?
557 +    //OOPSE = Object-Obfuscated Parallel Simulation Engine
558 +    
559 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, exclude_->getExcludeList(),
560 +                  &nGlobalExcludes, globalExcludes, molMembershipArray,
561 +                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
562 +
563 +    if( isError ){
564 +
565 +        sprintf( painCave.errMsg,
566 +                 "There was an error setting the simulation information in fortran.\n" );
567 +        painCave.isFatal = 1;
568 +        painCave.severity = OOPSE_ERROR;
569 +        simError();
570 +    }
571 +
572 + #ifdef IS_MPI
573 +    sprintf( checkPointMsg,
574 +       "succesfully sent the simulation information to fortran.\n");
575 +    MPIcheckPoint();
576 + #endif // is_mpi
577 + }
578 +
579 +
580 + #ifdef IS_MPI
581 + void SimInfo::setupFortranParallel() {
582 +    
583 +    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
584 +    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
585 +    std::vector<int> localToGlobalCutoffGroupIndex;
586 +    typename SimInfo::MoleculeIterator mi;
587 +    typename Molecule::AtomIterator ai;
588 +    typename Molecule::CutoffGroupIterator ci;
589 +    Molecule* mol;
590 +    Atom* atom;
591 +    CutoffGroup* cg;
592 +    mpiSimData parallelData;
593 +    int isError;
594 +
595 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
596 +
597 +        //local index(index in DataStorge) of atom is important
598 +        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
599 +            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
600 +        }
601 +
602 +        //local index of cutoff group is trivial, it only depends on the order of travesing
603 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
604 +            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
605 +        }        
606 +        
607 +    }
608 +
609 +    //fill up mpiSimData struct
610 +    parallelData.nMolGlobal = getNGlobalMolecules();
611 +    parallelData.nMolLocal = getNMolecules();
612 +    parallelData.nAtomsGlobal = getNGlobalAtoms();
613 +    parallelData.nAtomsLocal = getNAtoms();
614 +    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
615 +    parallelData.nGroupsLocal = getNCutoffGroups();
616 +    parallelData.myNode = worldRank;
617 +    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
618 +
619 +    //pass mpiSimData struct and index arrays to fortran
620 +    setFsimParallel(parallelData, &(parallelData->nAtomsLocal),
621 +                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
622 +                    &localToGlobalCutoffGroupIndex[0], &isError);
623 +
624 +    if (isError) {
625 +        sprintf(painCave.errMsg,
626 +                "mpiRefresh errror: fortran didn't like something we gave it.\n");
627 +        painCave.isFatal = 1;
628 +        simError();
629 +    }
630 +
631 +    sprintf(checkPointMsg, " mpiRefresh successful.\n");
632 +    MPIcheckPoint();
633 +
634 +
635 + }
636 +
637 + #endif
638 +
639 + double SimInfo::calcMaxCutoffRadius() {
640 +
641 +
642 +    std::vector<AtomType*> atomTypes;
643 +    std::vector<AtomType*>::iterator i;
644 +    std::vector<double> cutoffRadius;
645 +
646 +    //get the unique atom types
647 +    atomTypes = getUniqueAtomTypes();
648 +
649 +    //query the max cutoff radius among these atom types
650 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
651 +        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
652 +    }
653 +
654 +    double maxCutoffRadius = std::max_element(cutoffRadius.begin(), cutoffRadius.end());
655 + #ifdef IS_MPI
656 +    //pick the max cutoff radius among the processors
657 + #endif
658 +
659 +    return maxCutoffRadius;
660 + }
661 +
662 + void SimInfo::addProperty(GenericData* genData) {
663 +    properties_.addProperty(genData);  
664 + }
665 +
666 + void SimInfo::removeProperty(const std::string& propName) {
667 +    properties_.removeProperty(propName);  
668 + }
669 +
670 + void SimInfo::clearProperties() {
671 +    properties_.clearProperties();
672 + }
673 +
674 + std::vector<std::string> SimInfo::getPropertyNames() {
675 +    return properties_.getPropertyNames();  
676 + }
677 +      
678 + std::vector<GenericData*> SimInfo::getProperties() {
679 +    return properties_.getProperties();
680 + }
681 +
682 + GenericData* SimInfo::getPropertyByName(const std::string& propName) {
683 +    return properties_.getPropertyByName(propName);
684 + }
685 +
686 +
687 + std::ostream& operator <<(ostream& o, SimInfo& info) {
688 +
689 +    return o;
690 + }
691 +
692   }//end namespace oopse

Diff Legend

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