ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-2.0/src/brains/SimInfo.cpp
Revision: 2297
Committed: Thu Sep 15 00:14:35 2005 UTC (18 years, 10 months ago) by chrisfen
File size: 34847 byte(s)
Log Message:
changes to include the coulombicCorrection selector

File Contents

# User Rev Content
1 gezelter 2204 /*
2 gezelter 1930 * 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 gezelter 1490
49 gezelter 1930 #include <algorithm>
50     #include <set>
51 gezelter 1490
52 tim 1492 #include "brains/SimInfo.hpp"
53 gezelter 1930 #include "math/Vector3.hpp"
54     #include "primitives/Molecule.hpp"
55 gezelter 2285 #include "UseTheForce/fCutoffPolicy.h"
56 chrisfen 2297 #include "UseTheForce/fCoulombicCorrection.h"
57 gezelter 1930 #include "UseTheForce/doForces_interface.h"
58     #include "UseTheForce/notifyCutoffs_interface.h"
59     #include "utils/MemoryUtils.hpp"
60 tim 1492 #include "utils/simError.h"
61 tim 2000 #include "selection/SelectionManager.hpp"
62 gezelter 1490
63 gezelter 1930 #ifdef IS_MPI
64     #include "UseTheForce/mpiComponentPlan.h"
65     #include "UseTheForce/DarkSide/simParallel_interface.h"
66     #endif
67 gezelter 1490
68 gezelter 1930 namespace oopse {
69 gezelter 1490
70 gezelter 2204 SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs,
71     ForceField* ff, Globals* simParams) :
72     stamps_(stamps), forceField_(ff), simParams_(simParams),
73     ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
74     nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
75     nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
76     nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0),
77     nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0),
78     sman_(NULL), fortranInitialized_(false) {
79 gezelter 1490
80 gezelter 1930
81 gezelter 2204 std::vector<std::pair<MoleculeStamp*, int> >::iterator i;
82     MoleculeStamp* molStamp;
83     int nMolWithSameStamp;
84     int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
85     int nGroups = 0; //total cutoff groups defined in meta-data file
86     CutoffGroupStamp* cgStamp;
87     RigidBodyStamp* rbStamp;
88     int nRigidAtoms = 0;
89 gezelter 1930
90 gezelter 2204 for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) {
91 gezelter 1930 molStamp = i->first;
92     nMolWithSameStamp = i->second;
93    
94     addMoleculeStamp(molStamp, nMolWithSameStamp);
95 gezelter 1490
96 gezelter 1930 //calculate atoms in molecules
97     nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;
98 gezelter 1490
99    
100 gezelter 1930 //calculate atoms in cutoff groups
101     int nAtomsInGroups = 0;
102     int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
103    
104     for (int j=0; j < nCutoffGroupsInStamp; j++) {
105 gezelter 2204 cgStamp = molStamp->getCutoffGroup(j);
106     nAtomsInGroups += cgStamp->getNMembers();
107 gezelter 1930 }
108 gezelter 1490
109 gezelter 1930 nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
110     nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;
111 gezelter 1490
112 gezelter 1930 //calculate atoms in rigid bodies
113     int nAtomsInRigidBodies = 0;
114 tim 1958 int nRigidBodiesInStamp = molStamp->getNRigidBodies();
115 gezelter 1930
116     for (int j=0; j < nRigidBodiesInStamp; j++) {
117 gezelter 2204 rbStamp = molStamp->getRigidBody(j);
118     nAtomsInRigidBodies += rbStamp->getNMembers();
119 gezelter 1930 }
120 gezelter 1490
121 gezelter 1930 nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
122     nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;
123    
124 gezelter 2204 }
125 chrisfen 1636
126 gezelter 2204 //every free atom (atom does not belong to cutoff groups) is a cutoff group
127     //therefore the total number of cutoff groups in the system is equal to
128     //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data
129     //file plus the number of cutoff groups defined in meta-data file
130     nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
131 gezelter 1490
132 gezelter 2204 //every free atom (atom does not belong to rigid bodies) is an integrable object
133     //therefore the total number of integrable objects in the system is equal to
134     //the total number of atoms minus number of atoms belong to rigid body defined in meta-data
135     //file plus the number of rigid bodies defined in meta-data file
136     nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_;
137 gezelter 1490
138 gezelter 2204 nGlobalMols_ = molStampIds_.size();
139 gezelter 1490
140 gezelter 1930 #ifdef IS_MPI
141 gezelter 2204 molToProcMap_.resize(nGlobalMols_);
142 gezelter 1930 #endif
143 tim 1976
144 gezelter 2204 }
145 gezelter 1490
146 gezelter 2204 SimInfo::~SimInfo() {
147 tim 2082 std::map<int, Molecule*>::iterator i;
148     for (i = molecules_.begin(); i != molecules_.end(); ++i) {
149 gezelter 2204 delete i->second;
150 tim 2082 }
151     molecules_.clear();
152 tim 2187
153     delete stamps_;
154 gezelter 1930 delete sman_;
155     delete simParams_;
156     delete forceField_;
157 gezelter 2204 }
158 gezelter 1490
159 gezelter 2204 int SimInfo::getNGlobalConstraints() {
160 gezelter 1930 int nGlobalConstraints;
161     #ifdef IS_MPI
162     MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
163     MPI_COMM_WORLD);
164     #else
165     nGlobalConstraints = nConstraints_;
166     #endif
167     return nGlobalConstraints;
168 gezelter 2204 }
169 gezelter 1490
170 gezelter 2204 bool SimInfo::addMolecule(Molecule* mol) {
171 gezelter 1930 MoleculeIterator i;
172 gezelter 1490
173 gezelter 1930 i = molecules_.find(mol->getGlobalIndex());
174     if (i == molecules_.end() ) {
175 gezelter 1490
176 gezelter 2204 molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
177 gezelter 1930
178 gezelter 2204 nAtoms_ += mol->getNAtoms();
179     nBonds_ += mol->getNBonds();
180     nBends_ += mol->getNBends();
181     nTorsions_ += mol->getNTorsions();
182     nRigidBodies_ += mol->getNRigidBodies();
183     nIntegrableObjects_ += mol->getNIntegrableObjects();
184     nCutoffGroups_ += mol->getNCutoffGroups();
185     nConstraints_ += mol->getNConstraintPairs();
186 gezelter 1490
187 gezelter 2204 addExcludePairs(mol);
188 gezelter 1930
189 gezelter 2204 return true;
190 gezelter 1930 } else {
191 gezelter 2204 return false;
192 gezelter 1930 }
193 gezelter 2204 }
194 gezelter 1490
195 gezelter 2204 bool SimInfo::removeMolecule(Molecule* mol) {
196 gezelter 1930 MoleculeIterator i;
197     i = molecules_.find(mol->getGlobalIndex());
198 gezelter 1490
199 gezelter 1930 if (i != molecules_.end() ) {
200 gezelter 1490
201 gezelter 2204 assert(mol == i->second);
202 gezelter 1930
203 gezelter 2204 nAtoms_ -= mol->getNAtoms();
204     nBonds_ -= mol->getNBonds();
205     nBends_ -= mol->getNBends();
206     nTorsions_ -= mol->getNTorsions();
207     nRigidBodies_ -= mol->getNRigidBodies();
208     nIntegrableObjects_ -= mol->getNIntegrableObjects();
209     nCutoffGroups_ -= mol->getNCutoffGroups();
210     nConstraints_ -= mol->getNConstraintPairs();
211 gezelter 1490
212 gezelter 2204 removeExcludePairs(mol);
213     molecules_.erase(mol->getGlobalIndex());
214 gezelter 1490
215 gezelter 2204 delete mol;
216 gezelter 1930
217 gezelter 2204 return true;
218 gezelter 1930 } else {
219 gezelter 2204 return false;
220 gezelter 1930 }
221    
222    
223 gezelter 2204 }
224 gezelter 1930
225    
226 gezelter 2204 Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
227 gezelter 1930 i = molecules_.begin();
228     return i == molecules_.end() ? NULL : i->second;
229 gezelter 2204 }
230 gezelter 1930
231 gezelter 2204 Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
232 gezelter 1930 ++i;
233     return i == molecules_.end() ? NULL : i->second;
234 gezelter 2204 }
235 gezelter 1490
236    
237 gezelter 2204 void SimInfo::calcNdf() {
238 gezelter 1930 int ndf_local;
239     MoleculeIterator i;
240     std::vector<StuntDouble*>::iterator j;
241     Molecule* mol;
242     StuntDouble* integrableObject;
243 gezelter 1490
244 gezelter 1930 ndf_local = 0;
245    
246     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
247 gezelter 2204 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
248     integrableObject = mol->nextIntegrableObject(j)) {
249 gezelter 1490
250 gezelter 2204 ndf_local += 3;
251 gezelter 1490
252 gezelter 2204 if (integrableObject->isDirectional()) {
253     if (integrableObject->isLinear()) {
254     ndf_local += 2;
255     } else {
256     ndf_local += 3;
257     }
258     }
259 gezelter 1930
260 gezelter 2204 }//end for (integrableObject)
261 gezelter 1930 }// end for (mol)
262    
263     // n_constraints is local, so subtract them on each processor
264     ndf_local -= nConstraints_;
265    
266     #ifdef IS_MPI
267     MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
268     #else
269     ndf_ = ndf_local;
270     #endif
271    
272     // nZconstraints_ is global, as are the 3 COM translations for the
273     // entire system:
274     ndf_ = ndf_ - 3 - nZconstraint_;
275    
276 gezelter 2204 }
277 gezelter 1490
278 gezelter 2204 void SimInfo::calcNdfRaw() {
279 gezelter 1930 int ndfRaw_local;
280 gezelter 1490
281 gezelter 1930 MoleculeIterator i;
282     std::vector<StuntDouble*>::iterator j;
283     Molecule* mol;
284     StuntDouble* integrableObject;
285    
286     // Raw degrees of freedom that we have to set
287     ndfRaw_local = 0;
288    
289     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
290 gezelter 2204 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291     integrableObject = mol->nextIntegrableObject(j)) {
292 gezelter 1930
293 gezelter 2204 ndfRaw_local += 3;
294 gezelter 1930
295 gezelter 2204 if (integrableObject->isDirectional()) {
296     if (integrableObject->isLinear()) {
297     ndfRaw_local += 2;
298     } else {
299     ndfRaw_local += 3;
300     }
301     }
302 gezelter 1930
303 gezelter 2204 }
304 gezelter 1930 }
305    
306     #ifdef IS_MPI
307     MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
308     #else
309     ndfRaw_ = ndfRaw_local;
310     #endif
311 gezelter 2204 }
312 gezelter 1490
313 gezelter 2204 void SimInfo::calcNdfTrans() {
314 gezelter 1930 int ndfTrans_local;
315 gezelter 1490
316 gezelter 1930 ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
317 gezelter 1490
318    
319 gezelter 1930 #ifdef IS_MPI
320     MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
321     #else
322     ndfTrans_ = ndfTrans_local;
323     #endif
324 gezelter 1490
325 gezelter 1930 ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
326    
327 gezelter 2204 }
328 gezelter 1490
329 gezelter 2204 void SimInfo::addExcludePairs(Molecule* mol) {
330 gezelter 1930 std::vector<Bond*>::iterator bondIter;
331     std::vector<Bend*>::iterator bendIter;
332     std::vector<Torsion*>::iterator torsionIter;
333     Bond* bond;
334     Bend* bend;
335     Torsion* torsion;
336     int a;
337     int b;
338     int c;
339     int d;
340    
341     for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
342 gezelter 2204 a = bond->getAtomA()->getGlobalIndex();
343     b = bond->getAtomB()->getGlobalIndex();
344     exclude_.addPair(a, b);
345 gezelter 1930 }
346 gezelter 1490
347 gezelter 1930 for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
348 gezelter 2204 a = bend->getAtomA()->getGlobalIndex();
349     b = bend->getAtomB()->getGlobalIndex();
350     c = bend->getAtomC()->getGlobalIndex();
351 gezelter 1490
352 gezelter 2204 exclude_.addPair(a, b);
353     exclude_.addPair(a, c);
354     exclude_.addPair(b, c);
355 gezelter 1930 }
356 gezelter 1490
357 gezelter 1930 for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
358 gezelter 2204 a = torsion->getAtomA()->getGlobalIndex();
359     b = torsion->getAtomB()->getGlobalIndex();
360     c = torsion->getAtomC()->getGlobalIndex();
361     d = torsion->getAtomD()->getGlobalIndex();
362 gezelter 1490
363 gezelter 2204 exclude_.addPair(a, b);
364     exclude_.addPair(a, c);
365     exclude_.addPair(a, d);
366     exclude_.addPair(b, c);
367     exclude_.addPair(b, d);
368     exclude_.addPair(c, d);
369 gezelter 1490 }
370    
371 tim 2114 Molecule::RigidBodyIterator rbIter;
372     RigidBody* rb;
373     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
374 gezelter 2204 std::vector<Atom*> atoms = rb->getAtoms();
375     for (int i = 0; i < atoms.size() -1 ; ++i) {
376     for (int j = i + 1; j < atoms.size(); ++j) {
377     a = atoms[i]->getGlobalIndex();
378     b = atoms[j]->getGlobalIndex();
379     exclude_.addPair(a, b);
380     }
381     }
382 tim 2114 }
383    
384 gezelter 2204 }
385 gezelter 1930
386 gezelter 2204 void SimInfo::removeExcludePairs(Molecule* mol) {
387 gezelter 1930 std::vector<Bond*>::iterator bondIter;
388     std::vector<Bend*>::iterator bendIter;
389     std::vector<Torsion*>::iterator torsionIter;
390     Bond* bond;
391     Bend* bend;
392     Torsion* torsion;
393     int a;
394     int b;
395     int c;
396     int d;
397    
398     for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) {
399 gezelter 2204 a = bond->getAtomA()->getGlobalIndex();
400     b = bond->getAtomB()->getGlobalIndex();
401     exclude_.removePair(a, b);
402 gezelter 1490 }
403 gezelter 1930
404     for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) {
405 gezelter 2204 a = bend->getAtomA()->getGlobalIndex();
406     b = bend->getAtomB()->getGlobalIndex();
407     c = bend->getAtomC()->getGlobalIndex();
408 gezelter 1930
409 gezelter 2204 exclude_.removePair(a, b);
410     exclude_.removePair(a, c);
411     exclude_.removePair(b, c);
412 gezelter 1490 }
413 gezelter 1930
414     for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) {
415 gezelter 2204 a = torsion->getAtomA()->getGlobalIndex();
416     b = torsion->getAtomB()->getGlobalIndex();
417     c = torsion->getAtomC()->getGlobalIndex();
418     d = torsion->getAtomD()->getGlobalIndex();
419 gezelter 1930
420 gezelter 2204 exclude_.removePair(a, b);
421     exclude_.removePair(a, c);
422     exclude_.removePair(a, d);
423     exclude_.removePair(b, c);
424     exclude_.removePair(b, d);
425     exclude_.removePair(c, d);
426 gezelter 1930 }
427    
428 tim 2114 Molecule::RigidBodyIterator rbIter;
429     RigidBody* rb;
430     for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
431 gezelter 2204 std::vector<Atom*> atoms = rb->getAtoms();
432     for (int i = 0; i < atoms.size() -1 ; ++i) {
433     for (int j = i + 1; j < atoms.size(); ++j) {
434     a = atoms[i]->getGlobalIndex();
435     b = atoms[j]->getGlobalIndex();
436     exclude_.removePair(a, b);
437     }
438     }
439 tim 2114 }
440    
441 gezelter 2204 }
442 gezelter 1490
443    
444 gezelter 2204 void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
445 gezelter 1930 int curStampId;
446 gezelter 1490
447 gezelter 1930 //index from 0
448     curStampId = moleculeStamps_.size();
449 gezelter 1490
450 gezelter 1930 moleculeStamps_.push_back(molStamp);
451     molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
452 gezelter 2204 }
453 gezelter 1490
454 gezelter 2204 void SimInfo::update() {
455 gezelter 1490
456 gezelter 1930 setupSimType();
457 gezelter 1490
458 gezelter 1930 #ifdef IS_MPI
459     setupFortranParallel();
460     #endif
461 gezelter 1490
462 gezelter 1930 setupFortranSim();
463 gezelter 1490
464 gezelter 1930 //setup fortran force field
465     /** @deprecate */
466     int isError = 0;
467 chrisfen 2297
468     setupCoulombicCorrection( isError );
469    
470 gezelter 1930 if(isError){
471 gezelter 2204 sprintf( painCave.errMsg,
472     "ForceField error: There was an error initializing the forceField in fortran.\n" );
473     painCave.isFatal = 1;
474     simError();
475 gezelter 1930 }
476 gezelter 1490
477 gezelter 1930
478     setupCutoff();
479 gezelter 1490
480 gezelter 1930 calcNdf();
481     calcNdfRaw();
482     calcNdfTrans();
483    
484     fortranInitialized_ = true;
485 gezelter 2204 }
486 gezelter 1490
487 gezelter 2204 std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
488 gezelter 1930 SimInfo::MoleculeIterator mi;
489     Molecule* mol;
490     Molecule::AtomIterator ai;
491     Atom* atom;
492     std::set<AtomType*> atomTypes;
493 gezelter 1490
494 gezelter 1930 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
495 gezelter 1490
496 gezelter 2204 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
497     atomTypes.insert(atom->getAtomType());
498     }
499 gezelter 1930
500     }
501 gezelter 1490
502 gezelter 1930 return atomTypes;
503 gezelter 2204 }
504 gezelter 1490
505 gezelter 2204 void SimInfo::setupSimType() {
506 gezelter 1930 std::set<AtomType*>::iterator i;
507     std::set<AtomType*> atomTypes;
508     atomTypes = getUniqueAtomTypes();
509 gezelter 1490
510 gezelter 1930 int useLennardJones = 0;
511     int useElectrostatic = 0;
512     int useEAM = 0;
513     int useCharge = 0;
514     int useDirectional = 0;
515     int useDipole = 0;
516     int useGayBerne = 0;
517     int useSticky = 0;
518 chrisfen 2220 int useStickyPower = 0;
519 gezelter 1930 int useShape = 0;
520     int useFLARB = 0; //it is not in AtomType yet
521     int useDirectionalAtom = 0;
522     int useElectrostatics = 0;
523     //usePBC and useRF are from simParams
524     int usePBC = simParams_->getPBC();
525     int useRF = simParams_->getUseRF();
526 gezelter 1490
527 gezelter 1930 //loop over all of the atom types
528     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
529 gezelter 2204 useLennardJones |= (*i)->isLennardJones();
530     useElectrostatic |= (*i)->isElectrostatic();
531     useEAM |= (*i)->isEAM();
532     useCharge |= (*i)->isCharge();
533     useDirectional |= (*i)->isDirectional();
534     useDipole |= (*i)->isDipole();
535     useGayBerne |= (*i)->isGayBerne();
536     useSticky |= (*i)->isSticky();
537 chrisfen 2220 useStickyPower |= (*i)->isStickyPower();
538 gezelter 2204 useShape |= (*i)->isShape();
539 gezelter 1930 }
540 gezelter 1490
541 chrisfen 2220 if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
542 gezelter 2204 useDirectionalAtom = 1;
543 gezelter 1930 }
544 gezelter 1490
545 gezelter 1930 if (useCharge || useDipole) {
546 gezelter 2204 useElectrostatics = 1;
547 gezelter 1930 }
548 gezelter 1490
549 gezelter 1930 #ifdef IS_MPI
550     int temp;
551 gezelter 1490
552 gezelter 1930 temp = usePBC;
553     MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
554 gezelter 1490
555 gezelter 1930 temp = useDirectionalAtom;
556     MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
557 gezelter 1490
558 gezelter 1930 temp = useLennardJones;
559     MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
560 gezelter 1490
561 gezelter 1930 temp = useElectrostatics;
562     MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
563 gezelter 1490
564 gezelter 1930 temp = useCharge;
565     MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
566 gezelter 1490
567 gezelter 1930 temp = useDipole;
568     MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
569 gezelter 1490
570 gezelter 1930 temp = useSticky;
571     MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
572 gezelter 1490
573 chrisfen 2220 temp = useStickyPower;
574     MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
575    
576 gezelter 1930 temp = useGayBerne;
577     MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
578 gezelter 1490
579 gezelter 1930 temp = useEAM;
580     MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
581 gezelter 1490
582 gezelter 1930 temp = useShape;
583     MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
584    
585     temp = useFLARB;
586     MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
587    
588     temp = useRF;
589     MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
590 chrisfen 2279
591     temp = useUW;
592     MPI_Allreduce(&temp, &useUW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
593    
594     temp = useDW;
595     MPI_Allreduce(&temp, &useDW, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
596 gezelter 1930
597 gezelter 1490 #endif
598    
599 gezelter 1930 fInfo_.SIM_uses_PBC = usePBC;
600     fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
601     fInfo_.SIM_uses_LennardJones = useLennardJones;
602     fInfo_.SIM_uses_Electrostatics = useElectrostatics;
603     fInfo_.SIM_uses_Charges = useCharge;
604     fInfo_.SIM_uses_Dipoles = useDipole;
605     fInfo_.SIM_uses_Sticky = useSticky;
606 chrisfen 2220 fInfo_.SIM_uses_StickyPower = useStickyPower;
607 gezelter 1930 fInfo_.SIM_uses_GayBerne = useGayBerne;
608     fInfo_.SIM_uses_EAM = useEAM;
609     fInfo_.SIM_uses_Shapes = useShape;
610     fInfo_.SIM_uses_FLARB = useFLARB;
611     fInfo_.SIM_uses_RF = useRF;
612 gezelter 1490
613 gezelter 1930 if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) {
614 gezelter 1490
615 gezelter 2204 if (simParams_->haveDielectric()) {
616     fInfo_.dielect = simParams_->getDielectric();
617     } else {
618     sprintf(painCave.errMsg,
619     "SimSetup Error: No Dielectric constant was set.\n"
620     "\tYou are trying to use Reaction Field without"
621     "\tsetting a dielectric constant!\n");
622     painCave.isFatal = 1;
623     simError();
624     }
625 gezelter 1930
626     } else {
627 gezelter 2204 fInfo_.dielect = 0.0;
628 gezelter 1930 }
629    
630 gezelter 2204 }
631 gezelter 1490
632 gezelter 2204 void SimInfo::setupFortranSim() {
633 gezelter 1930 int isError;
634     int nExclude;
635     std::vector<int> fortranGlobalGroupMembership;
636    
637     nExclude = exclude_.getSize();
638     isError = 0;
639 gezelter 1490
640 gezelter 1930 //globalGroupMembership_ is filled by SimCreator
641     for (int i = 0; i < nGlobalAtoms_; i++) {
642 gezelter 2204 fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
643 gezelter 1930 }
644 gezelter 1490
645 gezelter 1930 //calculate mass ratio of cutoff group
646     std::vector<double> mfact;
647     SimInfo::MoleculeIterator mi;
648     Molecule* mol;
649     Molecule::CutoffGroupIterator ci;
650     CutoffGroup* cg;
651     Molecule::AtomIterator ai;
652     Atom* atom;
653     double totalMass;
654    
655     //to avoid memory reallocation, reserve enough space for mfact
656     mfact.reserve(getNCutoffGroups());
657 gezelter 1490
658 gezelter 1930 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
659 gezelter 2204 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
660 gezelter 1490
661 gezelter 2204 totalMass = cg->getMass();
662     for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
663     mfact.push_back(atom->getMass()/totalMass);
664     }
665 gezelter 1490
666 gezelter 2204 }
667 gezelter 1930 }
668 gezelter 1490
669 gezelter 1930 //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
670     std::vector<int> identArray;
671 gezelter 1490
672 gezelter 1930 //to avoid memory reallocation, reserve enough space identArray
673     identArray.reserve(getNAtoms());
674    
675     for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
676 gezelter 2204 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
677     identArray.push_back(atom->getIdent());
678     }
679 gezelter 1930 }
680 gezelter 1490
681 gezelter 1930 //fill molMembershipArray
682     //molMembershipArray is filled by SimCreator
683     std::vector<int> molMembershipArray(nGlobalAtoms_);
684     for (int i = 0; i < nGlobalAtoms_; i++) {
685 gezelter 2204 molMembershipArray[i] = globalMolMembership_[i] + 1;
686 gezelter 1930 }
687    
688     //setup fortran simulation
689     int nGlobalExcludes = 0;
690     int* globalExcludes = NULL;
691     int* excludeList = exclude_.getExcludeList();
692     setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList ,
693 gezelter 2204 &nGlobalExcludes, globalExcludes, &molMembershipArray[0],
694     &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError);
695 gezelter 1490
696 gezelter 1930 if( isError ){
697 gezelter 1490
698 gezelter 2204 sprintf( painCave.errMsg,
699     "There was an error setting the simulation information in fortran.\n" );
700     painCave.isFatal = 1;
701     painCave.severity = OOPSE_ERROR;
702     simError();
703 gezelter 1930 }
704    
705     #ifdef IS_MPI
706     sprintf( checkPointMsg,
707 gezelter 2204 "succesfully sent the simulation information to fortran.\n");
708 gezelter 1930 MPIcheckPoint();
709     #endif // is_mpi
710 gezelter 2204 }
711 gezelter 1490
712    
713 gezelter 1930 #ifdef IS_MPI
714 gezelter 2204 void SimInfo::setupFortranParallel() {
715 gezelter 1930
716     //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
717     std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
718     std::vector<int> localToGlobalCutoffGroupIndex;
719     SimInfo::MoleculeIterator mi;
720     Molecule::AtomIterator ai;
721     Molecule::CutoffGroupIterator ci;
722     Molecule* mol;
723     Atom* atom;
724     CutoffGroup* cg;
725     mpiSimData parallelData;
726     int isError;
727 gezelter 1490
728 gezelter 1930 for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
729 gezelter 1490
730 gezelter 2204 //local index(index in DataStorge) of atom is important
731     for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
732     localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
733     }
734 gezelter 1490
735 gezelter 2204 //local index of cutoff group is trivial, it only depends on the order of travesing
736     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
737     localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
738     }
739 gezelter 1930
740     }
741 gezelter 1490
742 gezelter 1930 //fill up mpiSimData struct
743     parallelData.nMolGlobal = getNGlobalMolecules();
744     parallelData.nMolLocal = getNMolecules();
745     parallelData.nAtomsGlobal = getNGlobalAtoms();
746     parallelData.nAtomsLocal = getNAtoms();
747     parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
748     parallelData.nGroupsLocal = getNCutoffGroups();
749     parallelData.myNode = worldRank;
750     MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
751 gezelter 1490
752 gezelter 1930 //pass mpiSimData struct and index arrays to fortran
753     setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
754     &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal),
755     &localToGlobalCutoffGroupIndex[0], &isError);
756 gezelter 1490
757 gezelter 1930 if (isError) {
758 gezelter 2204 sprintf(painCave.errMsg,
759     "mpiRefresh errror: fortran didn't like something we gave it.\n");
760     painCave.isFatal = 1;
761     simError();
762 gezelter 1930 }
763 gezelter 1490
764 gezelter 1930 sprintf(checkPointMsg, " mpiRefresh successful.\n");
765     MPIcheckPoint();
766 gezelter 1490
767    
768 gezelter 2204 }
769 chrisfen 1636
770 gezelter 1930 #endif
771 chrisfen 1636
772 gezelter 2204 double SimInfo::calcMaxCutoffRadius() {
773 chrisfen 1636
774    
775 gezelter 1930 std::set<AtomType*> atomTypes;
776     std::set<AtomType*>::iterator i;
777     std::vector<double> cutoffRadius;
778 gezelter 1490
779 gezelter 1930 //get the unique atom types
780     atomTypes = getUniqueAtomTypes();
781    
782     //query the max cutoff radius among these atom types
783     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
784 gezelter 2204 cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i));
785 gezelter 1930 }
786    
787     double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end()));
788 gezelter 1490 #ifdef IS_MPI
789 gezelter 1930 //pick the max cutoff radius among the processors
790 gezelter 1490 #endif
791    
792 gezelter 1930 return maxCutoffRadius;
793 gezelter 2204 }
794 gezelter 1930
795 gezelter 2204 void SimInfo::getCutoff(double& rcut, double& rsw) {
796 gezelter 1490
797 gezelter 1930 if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
798    
799 gezelter 2204 if (!simParams_->haveRcut()){
800     sprintf(painCave.errMsg,
801 gezelter 1930 "SimCreator Warning: No value was set for the cutoffRadius.\n"
802     "\tOOPSE will use a default value of 15.0 angstroms"
803     "\tfor the cutoffRadius.\n");
804 gezelter 2204 painCave.isFatal = 0;
805     simError();
806     rcut = 15.0;
807     } else{
808     rcut = simParams_->getRcut();
809     }
810 gezelter 1930
811 gezelter 2204 if (!simParams_->haveRsw()){
812     sprintf(painCave.errMsg,
813 gezelter 1930 "SimCreator Warning: No value was set for switchingRadius.\n"
814     "\tOOPSE will use a default value of\n"
815     "\t0.95 * cutoffRadius for the switchingRadius\n");
816 gezelter 2204 painCave.isFatal = 0;
817     simError();
818     rsw = 0.95 * rcut;
819     } else{
820     rsw = simParams_->getRsw();
821     }
822 gezelter 1930
823     } else {
824 gezelter 2204 // if charge, dipole or reaction field is not used and the cutofff radius is not specified in
825     //meta-data file, the maximum cutoff radius calculated from forcefiled will be used
826 gezelter 1930
827 gezelter 2204 if (simParams_->haveRcut()) {
828     rcut = simParams_->getRcut();
829     } else {
830     //set cutoff radius to the maximum cutoff radius based on atom types in the whole system
831     rcut = calcMaxCutoffRadius();
832     }
833 gezelter 1930
834 gezelter 2204 if (simParams_->haveRsw()) {
835     rsw = simParams_->getRsw();
836     } else {
837     rsw = rcut;
838     }
839 gezelter 1930
840     }
841 gezelter 2204 }
842 tim 2010
843 gezelter 2285 void SimInfo::setupCutoff() {
844 tim 2010 getCutoff(rcut_, rsw_);
845 gezelter 1930 double rnblist = rcut_ + 1; // skin of neighbor list
846    
847     //Pass these cutoff radius etc. to fortran. This function should be called once and only once
848 gezelter 2285
849     int cp = TRADITIONAL_CUTOFF_POLICY;
850     if (simParams_->haveCutoffPolicy()) {
851     std::string myPolicy = simParams_->getCutoffPolicy();
852     if (myPolicy == "MIX") {
853     cp = MIX_CUTOFF_POLICY;
854     } else {
855     if (myPolicy == "MAX") {
856     cp = MAX_CUTOFF_POLICY;
857     } else {
858     if (myPolicy == "TRADITIONAL") {
859     cp = TRADITIONAL_CUTOFF_POLICY;
860     } else {
861     // throw error
862     sprintf( painCave.errMsg,
863     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
864     painCave.isFatal = 1;
865     simError();
866     }
867     }
868     }
869     }
870     notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp);
871 gezelter 2204 }
872 gezelter 1490
873 chrisfen 2297 void SimInfo::setupCoulombicCorrection( int isError ) {
874    
875     int errorOut;
876     int cc = NONE;
877     double alphaVal;
878    
879     errorOut = isError;
880    
881     if (simParams_->haveCoulombicCorrection()) {
882     std::string myCorrection = simParams_->getCoulombicCorrection();
883     if (myCorrection == "NONE") {
884     cc = NONE;
885     } else {
886     if (myCorrection == "UNDAMPED_WOLF") {
887     cc = UNDAMPED_WOLF;
888     } else {
889     if (myCorrection == "WOLF") {
890     cc = WOLF;
891     if (!simParams_->haveDampingAlpha()) {
892     //throw error
893     sprintf( painCave.errMsg,
894     "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used for the Wolf Coulombic Correction.", simParams_->getDampingAlpha());
895     painCave.isFatal = 0;
896     simError();
897     }
898     alphaVal = simParams_->getDampingAlpha();
899     } else {
900     if (myCorrection == "REACTION_FIELD") {
901     cc = REACTION_FIELD;
902     } else {
903     // throw error
904     sprintf( painCave.errMsg,
905     "SimInfo error: Unknown coulombicCorrection. (Input file specified %s .)\n\tcoulombicCorrection must be one of: \"none\", \"undamped_wolf\", \"wolf\", or \"reaction_field\".", myCorrection.c_str() );
906     painCave.isFatal = 1;
907     simError();
908     }
909     }
910     }
911     }
912     }
913     initFortranFF( &fInfo_.SIM_uses_RF, &cc, &alphaVal, &errorOut );
914     }
915    
916 gezelter 2204 void SimInfo::addProperty(GenericData* genData) {
917 gezelter 1930 properties_.addProperty(genData);
918 gezelter 2204 }
919 gezelter 1490
920 gezelter 2204 void SimInfo::removeProperty(const std::string& propName) {
921 gezelter 1930 properties_.removeProperty(propName);
922 gezelter 2204 }
923 gezelter 1490
924 gezelter 2204 void SimInfo::clearProperties() {
925 gezelter 1930 properties_.clearProperties();
926 gezelter 2204 }
927 gezelter 1490
928 gezelter 2204 std::vector<std::string> SimInfo::getPropertyNames() {
929 gezelter 1930 return properties_.getPropertyNames();
930 gezelter 2204 }
931 gezelter 1930
932 gezelter 2204 std::vector<GenericData*> SimInfo::getProperties() {
933 gezelter 1930 return properties_.getProperties();
934 gezelter 2204 }
935 gezelter 1490
936 gezelter 2204 GenericData* SimInfo::getPropertyByName(const std::string& propName) {
937 gezelter 1930 return properties_.getPropertyByName(propName);
938 gezelter 2204 }
939 gezelter 1490
940 gezelter 2204 void SimInfo::setSnapshotManager(SnapshotManager* sman) {
941 tim 2116 if (sman_ == sman) {
942 gezelter 2204 return;
943 tim 2116 }
944     delete sman_;
945 gezelter 1930 sman_ = sman;
946 gezelter 1490
947 gezelter 1930 Molecule* mol;
948     RigidBody* rb;
949     Atom* atom;
950     SimInfo::MoleculeIterator mi;
951     Molecule::RigidBodyIterator rbIter;
952     Molecule::AtomIterator atomIter;;
953    
954     for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
955    
956 gezelter 2204 for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
957     atom->setSnapshotManager(sman_);
958     }
959 gezelter 1930
960 gezelter 2204 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
961     rb->setSnapshotManager(sman_);
962     }
963 gezelter 1930 }
964 gezelter 1490
965 gezelter 2204 }
966 gezelter 1490
967 gezelter 2204 Vector3d SimInfo::getComVel(){
968 gezelter 1930 SimInfo::MoleculeIterator i;
969     Molecule* mol;
970 gezelter 1490
971 gezelter 1930 Vector3d comVel(0.0);
972     double totalMass = 0.0;
973 gezelter 1490
974 gezelter 1930
975     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
976 gezelter 2204 double mass = mol->getMass();
977     totalMass += mass;
978     comVel += mass * mol->getComVel();
979 gezelter 1930 }
980 gezelter 1490
981 gezelter 1930 #ifdef IS_MPI
982     double tmpMass = totalMass;
983     Vector3d tmpComVel(comVel);
984     MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
985     MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
986     #endif
987    
988     comVel /= totalMass;
989    
990     return comVel;
991 gezelter 2204 }
992 gezelter 1490
993 gezelter 2204 Vector3d SimInfo::getCom(){
994 gezelter 1930 SimInfo::MoleculeIterator i;
995     Molecule* mol;
996 gezelter 1490
997 gezelter 1930 Vector3d com(0.0);
998     double totalMass = 0.0;
999    
1000     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1001 gezelter 2204 double mass = mol->getMass();
1002     totalMass += mass;
1003     com += mass * mol->getCom();
1004 gezelter 1930 }
1005 gezelter 1490
1006     #ifdef IS_MPI
1007 gezelter 1930 double tmpMass = totalMass;
1008     Vector3d tmpCom(com);
1009     MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1010     MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1011 gezelter 1490 #endif
1012    
1013 gezelter 1930 com /= totalMass;
1014 gezelter 1490
1015 gezelter 1930 return com;
1016 gezelter 1490
1017 gezelter 2204 }
1018 gezelter 1930
1019 gezelter 2204 std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1020 gezelter 1930
1021     return o;
1022 gezelter 2204 }
1023 chuckv 2252
1024    
1025     /*
1026     Returns center of mass and center of mass velocity in one function call.
1027     */
1028    
1029     void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1030     SimInfo::MoleculeIterator i;
1031     Molecule* mol;
1032    
1033    
1034     double totalMass = 0.0;
1035    
1036 gezelter 1930
1037 chuckv 2252 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1038     double mass = mol->getMass();
1039     totalMass += mass;
1040     com += mass * mol->getCom();
1041     comVel += mass * mol->getComVel();
1042     }
1043    
1044     #ifdef IS_MPI
1045     double tmpMass = totalMass;
1046     Vector3d tmpCom(com);
1047     Vector3d tmpComVel(comVel);
1048     MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1049     MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1050     MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1051     #endif
1052    
1053     com /= totalMass;
1054     comVel /= totalMass;
1055     }
1056    
1057     /*
1058     Return intertia tensor for entire system and angular momentum Vector.
1059 chuckv 2256
1060    
1061     [ Ixx -Ixy -Ixz ]
1062     J =| -Iyx Iyy -Iyz |
1063     [ -Izx -Iyz Izz ]
1064 chuckv 2252 */
1065    
1066     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1067    
1068    
1069     double xx = 0.0;
1070     double yy = 0.0;
1071     double zz = 0.0;
1072     double xy = 0.0;
1073     double xz = 0.0;
1074     double yz = 0.0;
1075     Vector3d com(0.0);
1076     Vector3d comVel(0.0);
1077    
1078     getComAll(com, comVel);
1079    
1080     SimInfo::MoleculeIterator i;
1081     Molecule* mol;
1082    
1083     Vector3d thisq(0.0);
1084     Vector3d thisv(0.0);
1085    
1086     double thisMass = 0.0;
1087    
1088    
1089    
1090    
1091     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1092    
1093     thisq = mol->getCom()-com;
1094     thisv = mol->getComVel()-comVel;
1095     thisMass = mol->getMass();
1096     // Compute moment of intertia coefficients.
1097     xx += thisq[0]*thisq[0]*thisMass;
1098     yy += thisq[1]*thisq[1]*thisMass;
1099     zz += thisq[2]*thisq[2]*thisMass;
1100    
1101     // compute products of intertia
1102     xy += thisq[0]*thisq[1]*thisMass;
1103     xz += thisq[0]*thisq[2]*thisMass;
1104     yz += thisq[1]*thisq[2]*thisMass;
1105    
1106     angularMomentum += cross( thisq, thisv ) * thisMass;
1107    
1108     }
1109    
1110    
1111     inertiaTensor(0,0) = yy + zz;
1112     inertiaTensor(0,1) = -xy;
1113     inertiaTensor(0,2) = -xz;
1114     inertiaTensor(1,0) = -xy;
1115 chuckv 2256 inertiaTensor(1,1) = xx + zz;
1116 chuckv 2252 inertiaTensor(1,2) = -yz;
1117     inertiaTensor(2,0) = -xz;
1118     inertiaTensor(2,1) = -yz;
1119     inertiaTensor(2,2) = xx + yy;
1120    
1121     #ifdef IS_MPI
1122     Mat3x3d tmpI(inertiaTensor);
1123     Vector3d tmpAngMom;
1124     MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1125     MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1126     #endif
1127    
1128     return;
1129     }
1130    
1131     //Returns the angular momentum of the system
1132     Vector3d SimInfo::getAngularMomentum(){
1133    
1134     Vector3d com(0.0);
1135     Vector3d comVel(0.0);
1136     Vector3d angularMomentum(0.0);
1137    
1138     getComAll(com,comVel);
1139    
1140     SimInfo::MoleculeIterator i;
1141     Molecule* mol;
1142    
1143 chuckv 2256 Vector3d thisr(0.0);
1144     Vector3d thisp(0.0);
1145 chuckv 2252
1146 chuckv 2256 double thisMass;
1147 chuckv 2252
1148     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1149 chuckv 2256 thisMass = mol->getMass();
1150     thisr = mol->getCom()-com;
1151     thisp = (mol->getComVel()-comVel)*thisMass;
1152 chuckv 2252
1153 chuckv 2256 angularMomentum += cross( thisr, thisp );
1154    
1155 chuckv 2252 }
1156    
1157     #ifdef IS_MPI
1158     Vector3d tmpAngMom;
1159     MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD);
1160     #endif
1161    
1162     return angularMomentum;
1163     }
1164    
1165    
1166 gezelter 1930 }//end namespace oopse
1167