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 1804 by tim, Tue Nov 30 19:58:25 2004 UTC

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

Diff Legend

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