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