ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/SimInfo.cpp
Revision: 2114
Committed: Thu Mar 10 23:56:42 2005 UTC (19 years, 4 months ago) by tim
File size: 30509 byte(s)
Log Message:
adding exclude pairs for rigidbody and cutoff group

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Acknowledgement of the program authors must be made in any
10 * publication of scientific results based in part on use of the
11 * program. An acceptable form of acknowledgement is citation of
12 * the article in which the program was described (Matthew
13 * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 * Parallel Simulation Engine for Molecular Dynamics,"
16 * J. Comput. Chem. 26, pp. 252-271 (2005))
17 *
18 * 2. Redistributions of source code must retain the above copyright
19 * notice, this list of conditions and the following disclaimer.
20 *
21 * 3. Redistributions in binary form must reproduce the above copyright
22 * notice, this list of conditions and the following disclaimer in the
23 * documentation and/or other materials provided with the
24 * distribution.
25 *
26 * This software is provided "AS IS," without a warranty of any
27 * kind. All express or implied conditions, representations and
28 * warranties, including any implied warranty of merchantability,
29 * fitness for a particular purpose or non-infringement, are hereby
30 * excluded. The University of Notre Dame and its licensors shall not
31 * be liable for any damages suffered by licensee as a result of
32 * using, modifying or distributing the software or its
33 * derivatives. In no event will the University of Notre Dame or its
34 * licensors be liable for any lost revenue, profit or data, or for
35 * direct, indirect, special, consequential, incidental or punitive
36 * damages, however caused and regardless of the theory of liability,
37 * arising out of the use of or inability to use software, even if the
38 * University of Notre Dame has been advised of the possibility of
39 * such damages.
40 */
41
42 /**
43 * @file SimInfo.cpp
44 * @author tlin
45 * @date 11/02/2004
46 * @version 1.0
47 */
48
49 #include <algorithm>
50 #include <set>
51
52 #include "brains/SimInfo.hpp"
53 #include "math/Vector3.hpp"
54 #include "primitives/Molecule.hpp"
55 #include "UseTheForce/doForces_interface.h"
56 #include "UseTheForce/notifyCutoffs_interface.h"
57 #include "utils/MemoryUtils.hpp"
58 #include "utils/simError.h"
59 #include "selection/SelectionManager.hpp"
60
61 #ifdef IS_MPI
62 #include "UseTheForce/mpiComponentPlan.h"
63 #include "UseTheForce/DarkSide/simParallel_interface.h"
64 #endif
65
66 namespace oopse {
67
68 SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
69 ForceField* ff, Globals* simParams) :
70 forceField_(ff), simParams_(simParams),
71 ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
72 nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
73 nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
74 nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
75 nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0),
76 sman_(NULL), fortranInitialized_(false) {
77
78
79 std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
80 MoleculeStamp* molStamp;
81 int nMolWithSameStamp;
82 int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
83 int nGroups = 0; //total cutoff groups defined in meta-data file
84 CutoffGroupStamp* cgStamp;
85 RigidBodyStamp* rbStamp;
86 int nRigidAtoms = 0;
87
88 for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
89 molStamp = i->first;
90 nMolWithSameStamp = i->second;
91
92 addMoleculeStamp(molStamp, nMolWithSameStamp);
93
94 //calculate atoms in molecules
95 nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;
96
97
98 //calculate atoms in cutoff groups
99 int nAtomsInGroups = 0;
100 int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
101
102 for (int j=0; j < nCutoffGroupsInStamp; j++) {
103 cgStamp = molStamp->getCutoffGroup(j);
104 nAtomsInGroups += cgStamp->getNMembers();
105 }
106
107 nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
108 nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;
109
110 //calculate atoms in rigid bodies
111 int nAtomsInRigidBodies = 0;
112 int nRigidBodiesInStamp = molStamp->getNRigidBodies();
113
114 for (int j=0; j < nRigidBodiesInStamp; j++) {
115 rbStamp = molStamp->getRigidBody(j);
116 nAtomsInRigidBodies += rbStamp->getNMembers();
117 }
118
119 nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
120 nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;
121
122 }
123
124 //every free atom (atom does not belong to cutoff groups) is a cutoff group
125 //therefore the total number of cutoff groups in the system is equal to
126 //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
127 //file plus the number of cutoff groups defined in meta-data file
128 nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
129
130 //every free atom (atom does not belong to rigid bodies) is an integrable object
131 //therefore the total number of integrable objects in the system is equal to
132 //the total number of atoms minus number of atoms belong to rigid body defined in meta-data
133 //file plus the number of rigid bodies defined in meta-data file
134 nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
135
136 nGlobalMols_ = molStampIds_.size();
137
138 #ifdef IS_MPI
139 molToProcMap_.resize(nGlobalMols_);
140 #endif
141
142 }
143
144 SimInfo::~SimInfo() {
145 std::map<int, Molecule*>::iterator i;
146 for (i = molecules_.begin(); i != molecules_.end(); ++i) {
147 delete i->second;
148 }
149 molecules_.clear();
150
151 MemoryUtils::deletePointers(moleculeStamps_);
152
153 delete sman_;
154 delete simParams_;
155 delete forceField_;
156 }
157
158 int SimInfo::getNGlobalConstraints() {
159 int nGlobalConstraints;
160 #ifdef IS_MPI
161 MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
162 MPI_COMM_WORLD);
163 #else
164 nGlobalConstraints = nConstraints_;
165 #endif
166 return nGlobalConstraints;
167 }
168
169 bool SimInfo::addMolecule(Molecule* mol) {
170 MoleculeIterator i;
171
172 i = molecules_.find(mol->getGlobalIndex());
173 if (i == molecules_.end() ) {
174
175 molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
176
177 nAtoms_ += mol->getNAtoms();
178 nBonds_ += mol->getNBonds();
179 nBends_ += mol->getNBends();
180 nTorsions_ += mol->getNTorsions();
181 nRigidBodies_ += mol->getNRigidBodies();
182 nIntegrableObjects_ += mol->getNIntegrableObjects();
183 nCutoffGroups_ += mol->getNCutoffGroups();
184 nConstraints_ += mol->getNConstraintPairs();
185
186 addExcludePairs(mol);
187
188 return true;
189 } else {
190 return false;
191 }
192 }
193
194 bool SimInfo::removeMolecule(Molecule* mol) {
195 MoleculeIterator i;
196 i = molecules_.find(mol->getGlobalIndex());
197
198 if (i != molecules_.end() ) {
199
200 assert(mol == i->second);
201
202 nAtoms_ -= mol->getNAtoms();
203 nBonds_ -= mol->getNBonds();
204 nBends_ -= mol->getNBends();
205 nTorsions_ -= mol->getNTorsions();
206 nRigidBodies_ -= mol->getNRigidBodies();
207 nIntegrableObjects_ -= mol->getNIntegrableObjects();
208 nCutoffGroups_ -= mol->getNCutoffGroups();
209 nConstraints_ -= mol->getNConstraintPairs();
210
211 removeExcludePairs(mol);
212 molecules_.erase(mol->getGlobalIndex());
213
214 delete mol;
215
216 return true;
217 } else {
218 return false;
219 }
220
221
222 }
223
224
225 Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
226 i = molecules_.begin();
227 return i == molecules_.end() ? NULL : i->second;
228 }
229
230 Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
231 ++i;
232 return i == molecules_.end() ? NULL : i->second;
233 }
234
235
236 void SimInfo::calcNdf() {
237 int ndf_local;
238 MoleculeIterator i;
239 std::vector<StuntDouble*>::iterator j;
240 Molecule* mol;
241 StuntDouble* integrableObject;
242
243 ndf_local = 0;
244
245 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
246 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
247 integrableObject = mol->nextIntegrableObject(j)) {
248
249 ndf_local += 3;
250
251 if (integrableObject->isDirectional()) {
252 if (integrableObject->isLinear()) {
253 ndf_local += 2;
254 } else {
255 ndf_local += 3;
256 }
257 }
258
259 }//end for (integrableObject)
260 }// end for (mol)
261
262 // n_constraints is local, so subtract them on each processor
263 ndf_local -= nConstraints_;
264
265 #ifdef IS_MPI
266 MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
267 #else
268 ndf_ = ndf_local;
269 #endif
270
271 // nZconstraints_ is global, as are the 3 COM translations for the
272 // entire system:
273 ndf_ = ndf_ - 3 - nZconstraint_;
274
275 }
276
277 void SimInfo::calcNdfRaw() {
278 int ndfRaw_local;
279
280 MoleculeIterator i;
281 std::vector<StuntDouble*>::iterator j;
282 Molecule* mol;
283 StuntDouble* integrableObject;
284
285 // Raw degrees of freedom that we have to set
286 ndfRaw_local = 0;
287
288 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
289 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
290 integrableObject = mol->nextIntegrableObject(j)) {
291
292 ndfRaw_local += 3;
293
294 if (integrableObject->isDirectional()) {
295 if (integrableObject->isLinear()) {
296 ndfRaw_local += 2;
297 } else {
298 ndfRaw_local += 3;
299 }
300 }
301
302 }
303 }
304
305 #ifdef IS_MPI
306 MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
307 #else
308 ndfRaw_ = ndfRaw_local;
309 #endif
310 }
311
312 void SimInfo::calcNdfTrans() {
313 int ndfTrans_local;
314
315 ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
316
317
318 #ifdef IS_MPI
319 MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
320 #else
321 ndfTrans_ = ndfTrans_local;
322 #endif
323
324 ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
325
326 }
327
328 void SimInfo::addExcludePairs(Molecule* mol) {
329 std::vector<Bond*>::iterator bondIter;
330 std::vector<Bend*>::iterator bendIter;
331 std::vector<Torsion*>::iterator torsionIter;
332 Bond* bond;
333 Bend* bend;
334 Torsion* torsion;
335 int a;
336 int b;
337 int c;
338 int d;
339
340 for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
341 a = bond->getAtomA()->getGlobalIndex();
342 b = bond->getAtomB()->getGlobalIndex();
343 exclude_.addPair(a, b);
344 }
345
346 for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
347 a = bend->getAtomA()->getGlobalIndex();
348 b = bend->getAtomB()->getGlobalIndex();
349 c = bend->getAtomC()->getGlobalIndex();
350
351 exclude_.addPair(a, b);
352 exclude_.addPair(a, c);
353 exclude_.addPair(b, c);
354 }
355
356 for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
357 a = torsion->getAtomA()->getGlobalIndex();
358 b = torsion->getAtomB()->getGlobalIndex();
359 c = torsion->getAtomC()->getGlobalIndex();
360 d = torsion->getAtomD()->getGlobalIndex();
361
362 exclude_.addPair(a, b);
363 exclude_.addPair(a, c);
364 exclude_.addPair(a, d);
365 exclude_.addPair(b, c);
366 exclude_.addPair(b, d);
367 exclude_.addPair(c, d);
368 }
369
370 Molecule::RigidBodyIterator rbIter;
371 RigidBody* rb;
372 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
373 std::vector<Atom*> atoms = rb->getAtoms();
374 for (int i = 0; i < atoms.size() -1 ; ++i) {
375 for (int j = i + 1; j < atoms.size(); ++j) {
376 a = atoms[i]->getGlobalIndex();
377 b = atoms[j]->getGlobalIndex();
378 exclude_.addPair(a, b);
379 }
380 }
381 }
382
383 Molecule::CutoffGroupIterator cgIter;
384 CutoffGroup* cg;
385 for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
386 std::vector<Atom*> atoms = cg->getAtoms();
387 for (int i = 0; i < atoms.size() -1 ; ++i) {
388 for (int j = i + 1; j < atoms.size(); ++j) {
389 a = atoms[i]->getGlobalIndex();
390 b = atoms[j]->getGlobalIndex();
391 exclude_.addPair(a, b);
392 }
393 }
394 }
395
396 }
397
398 void SimInfo::removeExcludePairs(Molecule* mol) {
399 std::vector<Bond*>::iterator bondIter;
400 std::vector<Bend*>::iterator bendIter;
401 std::vector<Torsion*>::iterator torsionIter;
402 Bond* bond;
403 Bend* bend;
404 Torsion* torsion;
405 int a;
406 int b;
407 int c;
408 int d;
409
410 for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
411 a = bond->getAtomA()->getGlobalIndex();
412 b = bond->getAtomB()->getGlobalIndex();
413 exclude_.removePair(a, b);
414 }
415
416 for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
417 a = bend->getAtomA()->getGlobalIndex();
418 b = bend->getAtomB()->getGlobalIndex();
419 c = bend->getAtomC()->getGlobalIndex();
420
421 exclude_.removePair(a, b);
422 exclude_.removePair(a, c);
423 exclude_.removePair(b, c);
424 }
425
426 for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
427 a = torsion->getAtomA()->getGlobalIndex();
428 b = torsion->getAtomB()->getGlobalIndex();
429 c = torsion->getAtomC()->getGlobalIndex();
430 d = torsion->getAtomD()->getGlobalIndex();
431
432 exclude_.removePair(a, b);
433 exclude_.removePair(a, c);
434 exclude_.removePair(a, d);
435 exclude_.removePair(b, c);
436 exclude_.removePair(b, d);
437 exclude_.removePair(c, d);
438 }
439
440 Molecule::RigidBodyIterator rbIter;
441 RigidBody* rb;
442 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
443 std::vector<Atom*> atoms = rb->getAtoms();
444 for (int i = 0; i < atoms.size() -1 ; ++i) {
445 for (int j = i + 1; j < atoms.size(); ++j) {
446 a = atoms[i]->getGlobalIndex();
447 b = atoms[j]->getGlobalIndex();
448 exclude_.removePair(a, b);
449 }
450 }
451 }
452
453 Molecule::CutoffGroupIterator cgIter;
454 CutoffGroup* cg;
455 for (cg = mol->beginCutoffGroup(cgIter); cg != NULL; cg = mol->nextCutoffGroup(cgIter)) {
456 std::vector<Atom*> atoms = cg->getAtoms();
457 for (int i = 0; i < atoms.size() -1 ; ++i) {
458 for (int j = i + 1; j < atoms.size(); ++j) {
459 a = atoms[i]->getGlobalIndex();
460 b = atoms[j]->getGlobalIndex();
461 exclude_.removePair(a, b);
462 }
463 }
464 }
465
466 }
467
468
469 void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
470 int curStampId;
471
472 //index from 0
473 curStampId = moleculeStamps_.size();
474
475 moleculeStamps_.push_back(molStamp);
476 molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
477 }
478
479 void SimInfo::update() {
480
481 setupSimType();
482
483 #ifdef IS_MPI
484 setupFortranParallel();
485 #endif
486
487 setupFortranSim();
488
489 //setup fortran force field
490 /** @deprecate */
491 int isError = 0;
492 initFortranFF( &fInfo_.SIM_uses_RF , &isError );
493 if(isError){
494 sprintf( painCave.errMsg,
495 "ForceField error: There was an error initializing the forceField in fortran.\n" );
496 painCave.isFatal = 1;
497 simError();
498 }
499
500
501 setupCutoff();
502
503 calcNdf();
504 calcNdfRaw();
505 calcNdfTrans();
506
507 fortranInitialized_ = true;
508 }
509
510 std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
511 SimInfo::MoleculeIterator mi;
512 Molecule* mol;
513 Molecule::AtomIterator ai;
514 Atom* atom;
515 std::set<AtomType*> atomTypes;
516
517 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
518
519 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
520 atomTypes.insert(atom->getAtomType());
521 }
522
523 }
524
525 return atomTypes;
526 }
527
528 void SimInfo::setupSimType() {
529 std::set<AtomType*>::iterator i;
530 std::set<AtomType*> atomTypes;
531 atomTypes = getUniqueAtomTypes();
532
533 int useLennardJones = 0;
534 int useElectrostatic = 0;
535 int useEAM = 0;
536 int useCharge = 0;
537 int useDirectional = 0;
538 int useDipole = 0;
539 int useGayBerne = 0;
540 int useSticky = 0;
541 int useShape = 0;
542 int useFLARB = 0; //it is not in AtomType yet
543 int useDirectionalAtom = 0;
544 int useElectrostatics = 0;
545 //usePBC and useRF are from simParams
546 int usePBC = simParams_->getPBC();
547 int useRF = simParams_->getUseRF();
548
549 //loop over all of the atom types
550 for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
551 useLennardJones |= (*i)->isLennardJones();
552 useElectrostatic |= (*i)->isElectrostatic();
553 useEAM |= (*i)->isEAM();
554 useCharge |= (*i)->isCharge();
555 useDirectional |= (*i)->isDirectional();
556 useDipole |= (*i)->isDipole();
557 useGayBerne |= (*i)->isGayBerne();
558 useSticky |= (*i)->isSticky();
559 useShape |= (*i)->isShape();
560 }
561
562 if (useSticky || useDipole || useGayBerne || useShape) {
563 useDirectionalAtom = 1;
564 }
565
566 if (useCharge || useDipole) {
567 useElectrostatics = 1;
568 }
569
570 #ifdef IS_MPI
571 int temp;
572
573 temp = usePBC;
574 MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
575
576 temp = useDirectionalAtom;
577 MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
578
579 temp = useLennardJones;
580 MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
581
582 temp = useElectrostatics;
583 MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
584
585 temp = useCharge;
586 MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
587
588 temp = useDipole;
589 MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
590
591 temp = useSticky;
592 MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
593
594 temp = useGayBerne;
595 MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
596
597 temp = useEAM;
598 MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
599
600 temp = useShape;
601 MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
602
603 temp = useFLARB;
604 MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
605
606 temp = useRF;
607 MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
608
609 #endif
610
611 fInfo_.SIM_uses_PBC = usePBC;
612 fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
613 fInfo_.SIM_uses_LennardJones = useLennardJones;
614 fInfo_.SIM_uses_Electrostatics = useElectrostatics;
615 fInfo_.SIM_uses_Charges = useCharge;
616 fInfo_.SIM_uses_Dipoles = useDipole;
617 fInfo_.SIM_uses_Sticky = useSticky;
618 fInfo_.SIM_uses_GayBerne = useGayBerne;
619 fInfo_.SIM_uses_EAM = useEAM;
620 fInfo_.SIM_uses_Shapes = useShape;
621 fInfo_.SIM_uses_FLARB = useFLARB;
622 fInfo_.SIM_uses_RF = useRF;
623
624 if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
625
626 if (simParams_->haveDielectric()) {
627 fInfo_.dielect = simParams_->getDielectric();
628 } else {
629 sprintf(painCave.errMsg,
630 "SimSetup Error: No Dielectric constant was set.\n"
631 "\tYou are trying to use Reaction Field without"
632 "\tsetting a dielectric constant!\n");
633 painCave.isFatal = 1;
634 simError();
635 }
636
637 } else {
638 fInfo_.dielect = 0.0;
639 }
640
641 }
642
643 void SimInfo::setupFortranSim() {
644 int isError;
645 int nExclude;
646 std::vector<int> fortranGlobalGroupMembership;
647
648 nExclude = exclude_.getSize();
649 isError = 0;
650
651 //globalGroupMembership_ is filled by SimCreator
652 for (int i = 0; i < nGlobalAtoms_; i++) {
653 fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
654 }
655
656 //calculate mass ratio of cutoff group
657 std::vector<double> mfact;
658 SimInfo::MoleculeIterator mi;
659 Molecule* mol;
660 Molecule::CutoffGroupIterator ci;
661 CutoffGroup* cg;
662 Molecule::AtomIterator ai;
663 Atom* atom;
664 double totalMass;
665
666 //to avoid memory reallocation, reserve enough space for mfact
667 mfact.reserve(getNCutoffGroups());
668
669 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
670 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
671
672 totalMass = cg->getMass();
673 for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
674 mfact.push_back(atom->getMass()/totalMass);
675 }
676
677 }
678 }
679
680 //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
681 std::vector<int> identArray;
682
683 //to avoid memory reallocation, reserve enough space identArray
684 identArray.reserve(getNAtoms());
685
686 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
687 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
688 identArray.push_back(atom->getIdent());
689 }
690 }
691
692 //fill molMembershipArray
693 //molMembershipArray is filled by SimCreator
694 std::vector<int> molMembershipArray(nGlobalAtoms_);
695 for (int i = 0; i < nGlobalAtoms_; i++) {
696 molMembershipArray[i] = globalMolMembership_[i] + 1;
697 }
698
699 //setup fortran simulation
700 //gloalExcludes and molMembershipArray should go away (They are never used)
701 //why the hell fortran need to know molecule?
702 //OOPSE = Object-Obfuscated Parallel Simulation Engine
703 int nGlobalExcludes = 0;
704 int* globalExcludes = NULL;
705 int* excludeList = exclude_.getExcludeList();
706 setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
707 &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
708 &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
709
710 if( isError ){
711
712 sprintf( painCave.errMsg,
713 "There was an error setting the simulation information in fortran.\n" );
714 painCave.isFatal = 1;
715 painCave.severity = OOPSE_ERROR;
716 simError();
717 }
718
719 #ifdef IS_MPI
720 sprintf( checkPointMsg,
721 "succesfully sent the simulation information to fortran.\n");
722 MPIcheckPoint();
723 #endif // is_mpi
724 }
725
726
727 #ifdef IS_MPI
728 void SimInfo::setupFortranParallel() {
729
730 //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
731 std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
732 std::vector<int> localToGlobalCutoffGroupIndex;
733 SimInfo::MoleculeIterator mi;
734 Molecule::AtomIterator ai;
735 Molecule::CutoffGroupIterator ci;
736 Molecule* mol;
737 Atom* atom;
738 CutoffGroup* cg;
739 mpiSimData parallelData;
740 int isError;
741
742 for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
743
744 //local index(index in DataStorge) of atom is important
745 for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
746 localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
747 }
748
749 //local index of cutoff group is trivial, it only depends on the order of travesing
750 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
751 localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
752 }
753
754 }
755
756 //fill up mpiSimData struct
757 parallelData.nMolGlobal = getNGlobalMolecules();
758 parallelData.nMolLocal = getNMolecules();
759 parallelData.nAtomsGlobal = getNGlobalAtoms();
760 parallelData.nAtomsLocal = getNAtoms();
761 parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
762 parallelData.nGroupsLocal = getNCutoffGroups();
763 parallelData.myNode = worldRank;
764 MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
765
766 //pass mpiSimData struct and index arrays to fortran
767 setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
768 &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal),
769 &localToGlobalCutoffGroupIndex[0], &isError);
770
771 if (isError) {
772 sprintf(painCave.errMsg,
773 "mpiRefresh errror: fortran didn't like something we gave it.\n");
774 painCave.isFatal = 1;
775 simError();
776 }
777
778 sprintf(checkPointMsg, " mpiRefresh successful.\n");
779 MPIcheckPoint();
780
781
782 }
783
784 #endif
785
786 double SimInfo::calcMaxCutoffRadius() {
787
788
789 std::set<AtomType*> atomTypes;
790 std::set<AtomType*>::iterator i;
791 std::vector<double> cutoffRadius;
792
793 //get the unique atom types
794 atomTypes = getUniqueAtomTypes();
795
796 //query the max cutoff radius among these atom types
797 for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
798 cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
799 }
800
801 double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
802 #ifdef IS_MPI
803 //pick the max cutoff radius among the processors
804 #endif
805
806 return maxCutoffRadius;
807 }
808
809 void SimInfo::getCutoff(double& rcut, double& rsw) {
810
811 if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
812
813 if (!simParams_->haveRcut()){
814 sprintf(painCave.errMsg,
815 "SimCreator Warning: No value was set for the cutoffRadius.\n"
816 "\tOOPSE will use a default value of 15.0 angstroms"
817 "\tfor the cutoffRadius.\n");
818 painCave.isFatal = 0;
819 simError();
820 rcut = 15.0;
821 } else{
822 rcut = simParams_->getRcut();
823 }
824
825 if (!simParams_->haveRsw()){
826 sprintf(painCave.errMsg,
827 "SimCreator Warning: No value was set for switchingRadius.\n"
828 "\tOOPSE will use a default value of\n"
829 "\t0.95 * cutoffRadius for the switchingRadius\n");
830 painCave.isFatal = 0;
831 simError();
832 rsw = 0.95 * rcut;
833 } else{
834 rsw = simParams_->getRsw();
835 }
836
837 } else {
838 // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
839 //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
840
841 if (simParams_->haveRcut()) {
842 rcut = simParams_->getRcut();
843 } else {
844 //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
845 rcut = calcMaxCutoffRadius();
846 }
847
848 if (simParams_->haveRsw()) {
849 rsw = simParams_->getRsw();
850 } else {
851 rsw = rcut;
852 }
853
854 }
855 }
856
857 void SimInfo::setupCutoff() {
858 getCutoff(rcut_, rsw_);
859 double rnblist = rcut_ + 1; // skin of neighbor list
860
861 //Pass these cutoff radius etc. to fortran. This function should be called once and only once
862 notifyFortranCutoffs(&rcut_, &rsw_, &rnblist);
863 }
864
865 void SimInfo::addProperty(GenericData* genData) {
866 properties_.addProperty(genData);
867 }
868
869 void SimInfo::removeProperty(const std::string& propName) {
870 properties_.removeProperty(propName);
871 }
872
873 void SimInfo::clearProperties() {
874 properties_.clearProperties();
875 }
876
877 std::vector<std::string> SimInfo::getPropertyNames() {
878 return properties_.getPropertyNames();
879 }
880
881 std::vector<GenericData*> SimInfo::getProperties() {
882 return properties_.getProperties();
883 }
884
885 GenericData* SimInfo::getPropertyByName(const std::string& propName) {
886 return properties_.getPropertyByName(propName);
887 }
888
889 void SimInfo::setSnapshotManager(SnapshotManager* sman) {
890 //if (sman_ == sman_) {
891 // return;
892 //}
893
894 //delete sman_;
895 sman_ = sman;
896
897 Molecule* mol;
898 RigidBody* rb;
899 Atom* atom;
900 SimInfo::MoleculeIterator mi;
901 Molecule::RigidBodyIterator rbIter;
902 Molecule::AtomIterator atomIter;;
903
904 for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
905
906 for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
907 atom->setSnapshotManager(sman_);
908 }
909
910 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
911 rb->setSnapshotManager(sman_);
912 }
913 }
914
915 }
916
917 Vector3d SimInfo::getComVel(){
918 SimInfo::MoleculeIterator i;
919 Molecule* mol;
920
921 Vector3d comVel(0.0);
922 double totalMass = 0.0;
923
924
925 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
926 double mass = mol->getMass();
927 totalMass += mass;
928 comVel += mass * mol->getComVel();
929 }
930
931 #ifdef IS_MPI
932 double tmpMass = totalMass;
933 Vector3d tmpComVel(comVel);
934 MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
935 MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
936 #endif
937
938 comVel /= totalMass;
939
940 return comVel;
941 }
942
943 Vector3d SimInfo::getCom(){
944 SimInfo::MoleculeIterator i;
945 Molecule* mol;
946
947 Vector3d com(0.0);
948 double totalMass = 0.0;
949
950 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
951 double mass = mol->getMass();
952 totalMass += mass;
953 com += mass * mol->getCom();
954 }
955
956 #ifdef IS_MPI
957 double tmpMass = totalMass;
958 Vector3d tmpCom(com);
959 MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
960 MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
961 #endif
962
963 com /= totalMass;
964
965 return com;
966
967 }
968
969 std::ostream& operator <<(std::ostream& o, SimInfo& info) {
970
971 return o;
972 }
973
974 }//end namespace oopse
975