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

Comparing branches/new_design/OOPSE-3.0/src/brains/SimInfo.cpp (file contents):
Revision 1710 by tim, Thu Nov 4 19:48:22 2004 UTC vs.
Revision 1824 by tim, Thu Dec 2 03:12:25 2004 UTC

# Line 31 | Line 31
31   */
32  
33   #include <algorithm>
34 <
34 > #include <set>
35   #include "brains/SimInfo.hpp"
36 + #include "primitives/Molecule.hpp"
37 + #include "UseTheForce/doForces_interface.h"
38 + #include "UseTheForce/notifyCutoffs_interface.h"
39   #include "utils/MemoryUtils.hpp"
40 + #include "utils/simError.h"
41  
42   namespace oopse {
43  
44 < SimInfo::SimInfo() : nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
45 <        nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), sman_(NULL){
44 > SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
45 >                                ForceField* ff, Globals* globals) :
46 >                                forceField_(ff), globals_(globals), nAtoms_(0), nBonds_(0),
47 >                                nBends_(0), nTorsions_(0), nRigidBodies_(0), nIntegrableObjects_(0),
48 >                                nCutoffGroups_(0), nConstraints_(0), nZconstraint_(0), sman_(NULL),
49 >                                fortranInitialized_(false) {
50  
51 +    std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
52 +    MoleculeStamp* molStamp;
53 +    int nMolWithSameStamp;
54 +    int nCutoffAtoms; // number of atoms belong to cutoff groups
55 +    int nGroups;          //total cutoff groups defined in meta-data file
56 +    CutoffGroupStamp* cgStamp;
57 +    int nAtomsInGroups;
58 +    int nCutoffGroupsInStamp;
59 +    
60 +    RigidBodyStamp* rbStamp;
61 +    int nAtomsInRigidBodies;
62 +    int nRigidBodiesInStamp;
63 +    int nRigidAtoms;
64 +    int nRigidBodies;
65 +        
66 +    nGlobalAtoms_ =  0;
67 +    
68 +    nGroups = 0;
69 +    nCutoffAtoms = 0;
70 +    nRigidBodies = 0;
71 +    
72 +    for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
73 +        molStamp = i->first;
74 +        nMolWithSameStamp = i->second;
75 +        
76 +        addMoleculeStamp(molStamp, nMolWithSameStamp);
77 +
78 +        //calculate atoms in molecules
79 +        nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;  
80 +
81 +
82 +        //calculate atoms in cutoff groups
83 +        nAtomsInGroups = 0;
84 +        nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
85 +        
86 +        for (int j=0; j < nCutoffGroupsInStamp; j++) {
87 +            cgStamp = molStamp->getCutoffGroup(j);
88 +            nAtomsInGroups += cgStamp->getNMembers();
89 +        }
90 +
91 +        nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
92 +        nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;            
93 +
94 +        //calculate atoms in rigid bodies
95 +        nAtomsInRigidBodies = 0;
96 +        nRigidBodiesInStamp = molStamp->getNCutoffGroups();
97 +        
98 +        for (int j=0; j < nRigidBodiesInStamp; j++) {
99 +            rbStamp = molStamp->getRigidBody(j);
100 +            nRigidBodiesInStamp += rbStamp->getNMembers();
101 +        }
102 +
103 +        nRigidBodies += nRigidBodiesInStamp * nMolWithSameStamp;
104 +        nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;            
105 +        
106 +    }
107 +
108 +    //every free atom (atom does not belong to cutoff groups) is a cutoff group
109 +    //therefore the total number of cutoff groups in the system is equal to
110 +    //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
111 +    //file plus the number of cutoff groups defined in meta-data file
112 +    nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
113 +
114 +    //every free atom (atom does not belong to rigid bodies) is a rigid body
115 +    //therefore the total number of cutoff groups in the system is equal to
116 +    //the total number of atoms minus number of atoms belong to  rigid body defined in meta-data
117 +    //file plus the number of  rigid bodies defined in meta-data file
118 +    nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nRigidBodies;
119 +
120 +    //initialize globalGroupMembership_, every element of this array will be 0
121 +    globalGroupMembership_.insert(globalGroupMembership_.end(), nGlobalAtoms_, 0);
122 +
123 +    nGlobalMols_ = molStampIds_.size();
124 +
125 + #ifdef IS_MPI    
126 +    molToProcMap_.resize(nGlobalMols_);
127 + #endif
128 +    
129   }
130  
131   SimInfo::~SimInfo() {
132 <    MemoryUtils::deleteVectorOfPointer(molecules_);
132 >    //MemoryUtils::deleteVectorOfPointer(molecules_);
133 >
134 >    MemoryUtils::deleteVectorOfPointer(moleculeStamps_);
135 >    
136      delete sman_;
137 +    delete globals_;
138 +    delete forceField_;
139  
140   }
141  
142  
143   bool SimInfo::addMolecule(Molecule* mol) {
144 <    std::vector<Molecule*>::iterator i;
145 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
144 >    MoleculeIterator i;
145 >
146 >    i = molecules_.find(mol->getGlobalIndex());
147      if (i != molecules_.end() ) {
56 –        molecules_.push_back(mol);
148  
149 +        molecules_.insert(make_pair(mol->getGlobalIndex(), mol));
150 +        
151          nAtoms_ += mol->getNAtoms();
152          nBonds_ += mol->getNBonds();
153          nBends_ += mol->getNBends();
# Line 71 | Line 164 | bool SimInfo::removeMolecule(Molecule* mol) {
164   }
165  
166   bool SimInfo::removeMolecule(Molecule* mol) {
167 <    std::vector<Molecule*>::iterator i;
168 <    i = std::find(molecules_.begin(), molecules_.end(), mol);
167 >    MoleculeIterator i;
168 >    i = molecules_.find(mol->getGlobalIndex());
169  
170      if (i != molecules_.end() ) {
171 <        molecules_.push_back(mol);
171 >
172 >        assert(mol == i->second);
173 >        
174          nAtoms_ -= mol->getNAtoms();
175          nBonds_ -= mol->getNBonds();
176          nBends_ -= mol->getNBends();
# Line 85 | Line 180 | bool SimInfo::removeMolecule(Molecule* mol) {
180          nCutoffGroups_ -= mol->getNCutoffGroups();
181          nConstraints_ -= mol->getNConstraints();
182  
183 +        molecules_.erase(mol->getGlobalIndex());
184 +
185 +        delete mol;
186 +        
187          return true;
188      } else {
189          return false;
# Line 94 | Line 193 | Molecule* SimInfo::beginMolecule(std::vector<Molecule*
193   }    
194  
195          
196 < Molecule* SimInfo::beginMolecule(std::vector<Molecule*>::iterator& i) {
196 > Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
197      i = molecules_.begin();
198 <    return i == molecules_.end() ? NULL : *i;
198 >    return i == molecules_.end() ? NULL : i->second;
199   }    
200  
201 < Molecule* SimInfo::nextMolecule(std::vector<Molecule*>::iterator& i) {
201 > Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
202      ++i;
203 <    return i == molecules_.end() ? NULL : *i;    
203 >    return i == molecules_.end() ? NULL : i->second;    
204   }
205  
206  
207 < void SimInfo::calcNDF() {
207 > void SimInfo::calcNdf() {
208      int ndf_local;
209 <    std::vector<Molecule*>::iterator i;
209 >    MoleculeIterator i;
210      std::vector<StuntDouble*>::iterator j;
211      Molecule* mol;
212      StuntDouble* integrableObject;
# Line 140 | Line 239 | void SimInfo::calcNDF() {
239      ndf_ = ndf_local;
240   #endif
241  
242 <    // nZconstraints is global, as are the 3 COM translations for the
242 >    // nZconstraints_ is global, as are the 3 COM translations for the
243      // entire system:
244 <    ndf_ = ndf_ - 3 - nZconstraints;
244 >    ndf_ = ndf_ - 3 - nZconstraint_;
245  
246   }
247  
248 < void SimInfo::calcNDFRaw() {
248 > void SimInfo::calcNdfRaw() {
249      int ndfRaw_local;
250  
251 <    std::vector<Molecule*>::iterator i;
251 >    MoleculeIterator i;
252      std::vector<StuntDouble*>::iterator j;
253      Molecule* mol;
254      StuntDouble* integrableObject;
# Line 181 | Line 280 | void SimInfo::calcNDFTrans() {
280   #endif
281   }
282  
283 < void SimInfo::calcNDFTrans() {
283 > void SimInfo::calcNdfTrans() {
284      int ndfTrans_local;
285  
286      ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
# Line 193 | Line 292 | void SimInfo::calcNDFTrans() {
292      ndfTrans_ = ndfTrans_local;
293   #endif
294  
295 <    ndfTrans_ = ndfTrans_ - 3 - nZconstraints;
295 >    ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
296  
297   }
298  
299 + void SimInfo::addExcludePairs(Molecule* mol) {
300 +    std::vector<Bond*>::iterator bondIter;
301 +    std::vector<Bend*>::iterator bendIter;
302 +    std::vector<Torsion*>::iterator torsionIter;
303 +    Bond* bond;
304 +    Bend* bend;
305 +    Torsion* torsion;
306 +    int a;
307 +    int b;
308 +    int c;
309 +    int d;
310 +    
311 +    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
312 +        a = bond->getAtomA()->getGlobalIndex();
313 +        b = bond->getAtomB()->getGlobalIndex();        
314 +        exclude_.addPair(a, b);
315 +    }
316  
317 +    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
318 +        a = bend->getAtomA()->getGlobalIndex();
319 +        b = bend->getAtomB()->getGlobalIndex();        
320 +        c = bend->getAtomC()->getGlobalIndex();
321 +
322 +        exclude_.addPair(a, b);
323 +        exclude_.addPair(a, c);
324 +        exclude_.addPair(b, c);        
325 +    }
326 +
327 +    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
328 +        a = torsion->getAtomA()->getGlobalIndex();
329 +        b = torsion->getAtomB()->getGlobalIndex();        
330 +        c = torsion->getAtomC()->getGlobalIndex();        
331 +        d = torsion->getAtomD()->getGlobalIndex();        
332 +
333 +        exclude_.addPair(a, b);
334 +        exclude_.addPair(a, c);
335 +        exclude_.addPair(a, d);
336 +        exclude_.addPair(b, c);
337 +        exclude_.addPair(b, d);
338 +        exclude_.addPair(c, d);        
339 +    }
340 +
341 +    
342 + }
343 +
344 + void SimInfo::removeExcludePairs(Molecule* mol) {
345 +    std::vector<Bond*>::iterator bondIter;
346 +    std::vector<Bend*>::iterator bendIter;
347 +    std::vector<Torsion*>::iterator torsionIter;
348 +    Bond* bond;
349 +    Bend* bend;
350 +    Torsion* torsion;
351 +    int a;
352 +    int b;
353 +    int c;
354 +    int d;
355 +    
356 +    for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
357 +        a = bond->getAtomA()->getGlobalIndex();
358 +        b = bond->getAtomB()->getGlobalIndex();        
359 +        exclude_.removePair(a, b);
360 +    }
361 +
362 +    for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
363 +        a = bend->getAtomA()->getGlobalIndex();
364 +        b = bend->getAtomB()->getGlobalIndex();        
365 +        c = bend->getAtomC()->getGlobalIndex();
366 +
367 +        exclude_.removePair(a, b);
368 +        exclude_.removePair(a, c);
369 +        exclude_.removePair(b, c);        
370 +    }
371 +
372 +    for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
373 +        a = torsion->getAtomA()->getGlobalIndex();
374 +        b = torsion->getAtomB()->getGlobalIndex();        
375 +        c = torsion->getAtomC()->getGlobalIndex();        
376 +        d = torsion->getAtomD()->getGlobalIndex();        
377 +
378 +        exclude_.removePair(a, b);
379 +        exclude_.removePair(a, c);
380 +        exclude_.removePair(a, d);
381 +        exclude_.removePair(b, c);
382 +        exclude_.removePair(b, d);
383 +        exclude_.removePair(c, d);        
384 +    }
385 +
386 + }
387 +
388 +
389 + void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
390 +    int curStampId;
391 +
392 +    //index from 0
393 +    curStampId = molStampIds_.size();
394 +
395 +    moleculeStamps_.push_back(molStamp);
396 +    molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
397 + }
398 +
399 + void SimInfo::update() {
400 +
401 +    setupSimType();
402 +
403 + #ifdef IS_MPI
404 +    setupFortranParallel();
405 + #endif
406 +
407 +    setupFortranSim();
408 +
409 +    //setup fortran force field
410 +    /** @deprecate */    
411 +    int isError = 0;
412 +    initFortranFF( &fInfo_.SIM_uses_RF , &isError );
413 +    if(isError){
414 +        sprintf( painCave.errMsg,
415 +         "ForceField error: There was an error initializing the forceField in fortran.\n" );
416 +        painCave.isFatal = 1;
417 +        simError();
418 +    }
419 +  
420 +    
421 +    setupCutoff();
422 +
423 +    calcNdf();
424 +    calcNdfRaw();
425 +    calcNdfTrans();
426 +
427 +    fortranInitialized_ = true;
428 + }
429 +
430 + std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
431 +    SimInfo::MoleculeIterator mi;
432 +    Molecule* mol;
433 +    Molecule::AtomIterator ai;
434 +    Atom* atom;
435 +    std::set<AtomType*> atomTypes;
436 +
437 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
438 +
439 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
440 +            atomTypes.insert(atom->getAtomType());
441 +        }
442 +        
443 +    }
444 +
445 +    return atomTypes;        
446 + }
447 +
448 + void SimInfo::setupSimType() {
449 +    std::set<AtomType*>::iterator i;
450 +    std::set<AtomType*> atomTypes;
451 +    atomTypes = getUniqueAtomTypes();
452 +    
453 +    int useLennardJones = 0;
454 +    int useElectrostatic = 0;
455 +    int useEAM = 0;
456 +    int useCharge = 0;
457 +    int useDirectional = 0;
458 +    int useDipole = 0;
459 +    int useGayBerne = 0;
460 +    int useSticky = 0;
461 +    int useShape = 0;
462 +    int useFLARB = 0; //it is not in AtomType yet
463 +    int useDirectionalAtom = 0;    
464 +    int useElectrostatics = 0;
465 +    //usePBC and useRF are from globals
466 +    bool usePBC = globals_->getPBC();
467 +    bool useRF = globals_->getUseRF();
468 +
469 +    //loop over all of the atom types
470 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
471 +        useLennardJones |= (*i)->isLennardJones();
472 +        useElectrostatic |= (*i)->isElectrostatic();
473 +        useEAM |= (*i)->isEAM();
474 +        useCharge |= (*i)->isCharge();
475 +        useDirectional |= (*i)->isDirectional();
476 +        useDipole |= (*i)->isDipole();
477 +        useGayBerne |= (*i)->isGayBerne();
478 +        useSticky |= (*i)->isSticky();
479 +        useShape |= (*i)->isShape();
480 +    }
481 +
482 +    if (useSticky || useDipole || useGayBerne || useShape) {
483 +        useDirectionalAtom = 1;
484 +    }
485 +
486 +    if (useCharge || useDipole) {
487 +        useElectrostatics = 1;
488 +    }
489 +
490 + #ifdef IS_MPI    
491 +    int temp;
492 +
493 +    temp = usePBC;
494 +    MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
495 +
496 +    temp = useDirectionalAtom;
497 +    MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
498 +
499 +    temp = useLennardJones;
500 +    MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
501 +
502 +    temp = useElectrostatics;
503 +    MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
504 +
505 +    temp = useCharge;
506 +    MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
507 +
508 +    temp = useDipole;
509 +    MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
510 +
511 +    temp = useSticky;
512 +    MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
513 +
514 +    temp = useGayBerne;
515 +    MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
516 +
517 +    temp = useEAM;
518 +    MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
519 +
520 +    temp = useShape;
521 +    MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);  
522 +
523 +    temp = useFLARB;
524 +    MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
525 +
526 +    temp = useRF;
527 +    MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);    
528 +    
529 + #endif
530 +
531 +    fInfo_.SIM_uses_PBC = usePBC;    
532 +    fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
533 +    fInfo_.SIM_uses_LennardJones = useLennardJones;
534 +    fInfo_.SIM_uses_Electrostatics = useElectrostatics;    
535 +    fInfo_.SIM_uses_Charges = useCharge;
536 +    fInfo_.SIM_uses_Dipoles = useDipole;
537 +    fInfo_.SIM_uses_Sticky = useSticky;
538 +    fInfo_.SIM_uses_GayBerne = useGayBerne;
539 +    fInfo_.SIM_uses_EAM = useEAM;
540 +    fInfo_.SIM_uses_Shapes = useShape;
541 +    fInfo_.SIM_uses_FLARB = useFLARB;
542 +    fInfo_.SIM_uses_RF = useRF;
543 +
544 +    if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
545 +
546 +        if (globals_->haveDielectric()) {
547 +            fInfo_.dielect = globals_->getDielectric();
548 +        } else {
549 +            sprintf(painCave.errMsg,
550 +                    "SimSetup Error: No Dielectric constant was set.\n"
551 +                    "\tYou are trying to use Reaction Field without"
552 +                    "\tsetting a dielectric constant!\n");
553 +            painCave.isFatal = 1;
554 +            simError();
555 +        }
556 +        
557 +    } else {
558 +        fInfo_.dielect = 0.0;
559 +    }
560 +
561 + }
562 +
563 + void SimInfo::setupFortranSim() {
564 +    int isError;
565 +    int nExclude;
566 +    std::vector<int> fortranGlobalGroupMembership;
567 +    
568 +    nExclude = exclude_.getSize();
569 +    isError = 0;
570 +
571 +    //globalGroupMembership_ is filled by SimCreator    
572 +    for (int i = 0; i < nGlobalAtoms_; i++) {
573 +        fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
574 +    }
575 +
576 +    //calculate mass ratio of cutoff group
577 +    std::vector<double> mfact;
578 +    SimInfo::MoleculeIterator mi;
579 +    Molecule* mol;
580 +    Molecule::CutoffGroupIterator ci;
581 +    CutoffGroup* cg;
582 +    Molecule::AtomIterator ai;
583 +    Atom* atom;
584 +    double totalMass;
585 +
586 +    //to avoid memory reallocation, reserve enough space for mfact
587 +    mfact.reserve(getNCutoffGroups());
588 +    
589 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
590 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
591 +
592 +            totalMass = cg->getMass();
593 +            for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
594 +                        mfact.push_back(atom->getMass()/totalMass);
595 +            }
596 +
597 +        }      
598 +    }
599 +
600 +    //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
601 +    std::vector<int> identArray;
602 +
603 +    //to avoid memory reallocation, reserve enough space identArray
604 +    identArray.reserve(getNAtoms());
605 +    
606 +    for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {        
607 +        for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
608 +            identArray.push_back(atom->getIdent());
609 +        }
610 +    }    
611 +
612 +    //fill molMembershipArray
613 +    //molMembershipArray is filled by SimCreator    
614 +    std::vector<int> molMembershipArray(nGlobalAtoms_);
615 +    for (int i = 0; i < nGlobalAtoms_; i++) {
616 +        molMembershipArray.push_back(globalMolMembership_[i] + 1);
617 +    }
618 +    
619 +    //setup fortran simulation
620 +    //gloalExcludes and molMembershipArray should go away (They are never used)
621 +    //why the hell fortran need to know molecule?
622 +    //OOPSE = Object-Obfuscated Parallel Simulation Engine
623 +    int nGlobalExcludes = 0;
624 +    int* globalExcludes = NULL;
625 +    int* excludeList = exclude_.getExcludeList();
626 +    setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
627 +                  &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
628 +                  &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
629 +
630 +    if( isError ){
631 +
632 +        sprintf( painCave.errMsg,
633 +                 "There was an error setting the simulation information in fortran.\n" );
634 +        painCave.isFatal = 1;
635 +        painCave.severity = OOPSE_ERROR;
636 +        simError();
637 +    }
638 +
639 + #ifdef IS_MPI
640 +    sprintf( checkPointMsg,
641 +       "succesfully sent the simulation information to fortran.\n");
642 +    MPIcheckPoint();
643 + #endif // is_mpi
644 + }
645 +
646 +
647 + #ifdef IS_MPI
648 + void SimInfo::setupFortranParallel() {
649 +    
650 +    //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
651 +    std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
652 +    std::vector<int> localToGlobalCutoffGroupIndex;
653 +    SimInfo::MoleculeIterator mi;
654 +    Molecule::AtomIterator ai;
655 +    Molecule::CutoffGroupIterator ci;
656 +    Molecule* mol;
657 +    Atom* atom;
658 +    CutoffGroup* cg;
659 +    mpiSimData parallelData;
660 +    int isError;
661 +
662 +    for (mol = beginMolecule(mi); mol != NULL; mol  = nextMolecule(mi)) {
663 +
664 +        //local index(index in DataStorge) of atom is important
665 +        for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
666 +            localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
667 +        }
668 +
669 +        //local index of cutoff group is trivial, it only depends on the order of travesing
670 +        for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
671 +            localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
672 +        }        
673 +        
674 +    }
675 +
676 +    //fill up mpiSimData struct
677 +    parallelData.nMolGlobal = getNGlobalMolecules();
678 +    parallelData.nMolLocal = getNMolecules();
679 +    parallelData.nAtomsGlobal = getNGlobalAtoms();
680 +    parallelData.nAtomsLocal = getNAtoms();
681 +    parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
682 +    parallelData.nGroupsLocal = getNCutoffGroups();
683 +    parallelData.myNode = worldRank;
684 +    MPI_Comm_size(MPI_COMM_WORLD, &(parallelData->nProcessors));
685 +
686 +    //pass mpiSimData struct and index arrays to fortran
687 +    setFsimParallel(parallelData, &(parallelData->nAtomsLocal),
688 +                    &localToGlobalAtomIndex[0],  &(parallelData->nGroupsLocal),
689 +                    &localToGlobalCutoffGroupIndex[0], &isError);
690 +
691 +    if (isError) {
692 +        sprintf(painCave.errMsg,
693 +                "mpiRefresh errror: fortran didn't like something we gave it.\n");
694 +        painCave.isFatal = 1;
695 +        simError();
696 +    }
697 +
698 +    sprintf(checkPointMsg, " mpiRefresh successful.\n");
699 +    MPIcheckPoint();
700 +
701 +
702 + }
703 +
704 + #endif
705 +
706 + double SimInfo::calcMaxCutoffRadius() {
707 +
708 +
709 +    std::set<AtomType*> atomTypes;
710 +    std::set<AtomType*>::iterator i;
711 +    std::vector<double> cutoffRadius;
712 +
713 +    //get the unique atom types
714 +    atomTypes = getUniqueAtomTypes();
715 +
716 +    //query the max cutoff radius among these atom types
717 +    for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
718 +        cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
719 +    }
720 +
721 +    double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
722 + #ifdef IS_MPI
723 +    //pick the max cutoff radius among the processors
724 + #endif
725 +
726 +    return maxCutoffRadius;
727 + }
728 +
729 + void SimInfo::setupCutoff() {
730 +    double rcut_;  //cutoff radius
731 +    double rsw_; //switching radius
732 +    
733 +    if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
734 +        
735 +        if (!globals_->haveRcut()){
736 +            sprintf(painCave.errMsg,
737 +                "SimCreator Warning: No value was set for the cutoffRadius.\n"
738 +                "\tOOPSE will use a default value of 15.0 angstroms"
739 +                "\tfor the cutoffRadius.\n");
740 +            painCave.isFatal = 0;
741 +            simError();
742 +            rcut_ = 15.0;
743 +        } else{
744 +            rcut_ = globals_->getRcut();
745 +        }
746 +
747 +        if (!globals_->haveRsw()){
748 +            sprintf(painCave.errMsg,
749 +                "SimCreator Warning: No value was set for switchingRadius.\n"
750 +                "\tOOPSE will use a default value of\n"
751 +                "\t0.95 * cutoffRadius for the switchingRadius\n");
752 +            painCave.isFatal = 0;
753 +            simError();
754 +            rsw_ = 0.95 * rcut_;
755 +        } else{
756 +            rsw_ = globals_->getRsw();
757 +        }
758 +
759 +    } else {
760 +        // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
761 +        //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
762 +        
763 +        if (globals_->haveRcut()) {
764 +            rcut_ = globals_->getRcut();
765 +        } else {
766 +            //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
767 +            rcut_ = calcMaxCutoffRadius();
768 +        }
769 +
770 +        if (globals_->haveRsw()) {
771 +            rsw_  = globals_->getRsw();
772 +        } else {
773 +            rsw_ = rcut_;
774 +        }
775 +    
776 +    }
777 +        
778 +    double rnblist = rcut_ + 1; // skin of neighbor list
779 +
780 +    //Pass these cutoff radius etc. to fortran. This function should be called once and only once
781 +    notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
782 + }
783 +
784 + void SimInfo::addProperty(GenericData* genData) {
785 +    properties_.addProperty(genData);  
786 + }
787 +
788 + void SimInfo::removeProperty(const std::string& propName) {
789 +    properties_.removeProperty(propName);  
790 + }
791 +
792 + void SimInfo::clearProperties() {
793 +    properties_.clearProperties();
794 + }
795 +
796 + std::vector<std::string> SimInfo::getPropertyNames() {
797 +    return properties_.getPropertyNames();  
798 + }
799 +      
800 + std::vector<GenericData*> SimInfo::getProperties() {
801 +    return properties_.getProperties();
802 + }
803 +
804 + GenericData* SimInfo::getPropertyByName(const std::string& propName) {
805 +    return properties_.getPropertyByName(propName);
806 + }
807 +
808 + void SimInfo::setSnapshotManager(SnapshotManager* sman) {
809 +    sman_ = sman;
810 +
811 +    Molecule* mol;
812 +    RigidBody* rb;
813 +    Atom* atom;
814 +    SimInfo::MoleculeIterator mi;
815 +    Molecule::RigidBodyIterator rbIter;
816 +    Molecule::AtomIterator atomIter;;
817 +
818 +    for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
819 +        
820 +        for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
821 +            atom->setSnapshotManager(sman_);
822 +        }
823 +        
824 +        for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
825 +            rb->setSnapshotManager(sman_);
826 +        }
827 +    }    
828 +    
829 + }
830 +
831 + std::ostream& operator <<(ostream& o, SimInfo& info) {
832 +
833 +    return o;
834 + }
835 +
836   }//end namespace oopse

Diff Legend

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