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_); |
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 |
|
MoleculeIterator i; |
109 |
< |
i = std::find(molecules_.begin(), molecules_.end(), mol); |
109 |
> |
|
110 |
> |
i = molecules_.find(mol->getGlobalIndex()); |
111 |
|
if (i != molecules_.end() ) { |
112 |
|
|
113 |
|
molecules_.insert(make_pair(mol->getGlobalIndex(), mol)); |
129 |
|
|
130 |
|
bool SimInfo::removeMolecule(Molecule* mol) { |
131 |
|
MoleculeIterator i; |
132 |
< |
i = std::find(molecules_.begin(), molecules_.end(), mol); |
132 |
> |
i = molecules_.find(mol->getGlobalIndex()); |
133 |
|
|
134 |
|
if (i != molecules_.end() ) { |
135 |
|
|
136 |
+ |
assert(mol == i->second); |
137 |
+ |
|
138 |
|
nAtoms_ -= mol->getNAtoms(); |
139 |
|
nBonds_ -= mol->getNBonds(); |
140 |
|
nBends_ -= mol->getNBends(); |
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 |
|
|
256 |
|
ndfTrans_ = ndfTrans_local; |
257 |
|
#endif |
258 |
|
|
259 |
< |
ndfTrans_ = ndfTrans_ - 3 - nZconstraints; |
259 |
> |
ndfTrans_ = ndfTrans_ - 3 - nZconstraints_; |
260 |
|
|
261 |
|
} |
262 |
|
|
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; |
606 |
|
|
607 |
|
} |
608 |
|
|
609 |
< |
//Setup Parallel Data and pass the index arrays to fortran |
610 |
< |
parallelData.nMolGlobal = getNMolGlobal(); |
611 |
< |
parallelData.nMolLocal = ; |
612 |
< |
parallelData.nAtomsGlobal = ; |
613 |
< |
parallelData.nAtomsLocal = ; |
614 |
< |
parallelData.nGroupsGlobal = ; |
615 |
< |
parallelData.nGroupsLocal = ; |
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 |
< |
setFsimParallel(parallelData, &(parallelData->nAtomsLocal), |
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 |
|
|
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; |