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

Diff Legend

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