ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/brains/SimInfo.cpp
Revision: 3559
Committed: Fri Oct 23 18:41:09 2009 UTC (15 years, 6 months ago) by gezelter
File size: 52046 byte(s)
Log Message:
removing MPI responsibilities from the lowest level force routines.  This is
in preparation for migrating these routines (LJ, electrostatic, eam, suttonchen,
gay-berne, sticky) to C++

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 tim 2448 #include <map>
52 gezelter 1490
53 tim 1492 #include "brains/SimInfo.hpp"
54 gezelter 1930 #include "math/Vector3.hpp"
55     #include "primitives/Molecule.hpp"
56 tim 2982 #include "primitives/StuntDouble.hpp"
57 gezelter 2285 #include "UseTheForce/fCutoffPolicy.h"
58 chrisfen 2305 #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h"
59 chrisfen 2415 #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h"
60 chrisfen 2425 #include "UseTheForce/DarkSide/fSwitchingFunctionType.h"
61 gezelter 1930 #include "UseTheForce/doForces_interface.h"
62 chuckv 3080 #include "UseTheForce/DarkSide/neighborLists_interface.h"
63 chrisfen 2309 #include "UseTheForce/DarkSide/electrostatic_interface.h"
64 chrisfen 2425 #include "UseTheForce/DarkSide/switcheroo_interface.h"
65 gezelter 1930 #include "utils/MemoryUtils.hpp"
66 tim 1492 #include "utils/simError.h"
67 tim 2000 #include "selection/SelectionManager.hpp"
68 chuckv 2533 #include "io/ForceFieldOptions.hpp"
69     #include "UseTheForce/ForceField.hpp"
70 gezelter 1490
71 chuckv 3080
72 gezelter 1930 #ifdef IS_MPI
73     #include "UseTheForce/mpiComponentPlan.h"
74     #include "UseTheForce/DarkSide/simParallel_interface.h"
75     #endif
76 gezelter 1490
77 gezelter 1930 namespace oopse {
78 tim 2448 std::set<int> getRigidSet(int index, std::map<int, std::set<int> >& container) {
79     std::map<int, std::set<int> >::iterator i = container.find(index);
80     std::set<int> result;
81     if (i != container.end()) {
82     result = i->second;
83     }
84 gezelter 1490
85 tim 2448 return result;
86     }
87    
88 tim 2469 SimInfo::SimInfo(ForceField* ff, Globals* simParams) :
89     forceField_(ff), simParams_(simParams),
90 gezelter 2733 ndf_(0), fdf_local(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0),
91 gezelter 2204 nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0),
92     nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0),
93 gezelter 3432 nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nInversions_(0),
94     nRigidBodies_(0), nIntegrableObjects_(0), nCutoffGroups_(0),
95     nConstraints_(0), sman_(NULL), fortranInitialized_(false),
96     calcBoxDipole_(false), useAtomicVirial_(true) {
97 gezelter 1490
98 gezelter 3432
99 gezelter 2204 MoleculeStamp* molStamp;
100     int nMolWithSameStamp;
101     int nCutoffAtoms = 0; // number of atoms belong to cutoff groups
102 chrisfen 2344 int nGroups = 0; //total cutoff groups defined in meta-data file
103 gezelter 2204 CutoffGroupStamp* cgStamp;
104     RigidBodyStamp* rbStamp;
105     int nRigidAtoms = 0;
106 gezelter 3432
107 tim 2469 std::vector<Component*> components = simParams->getComponents();
108    
109     for (std::vector<Component*>::iterator i = components.begin(); i !=components.end(); ++i) {
110     molStamp = (*i)->getMoleculeStamp();
111     nMolWithSameStamp = (*i)->getNMol();
112 gezelter 1930
113     addMoleculeStamp(molStamp, nMolWithSameStamp);
114 gezelter 1490
115 gezelter 1930 //calculate atoms in molecules
116     nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp;
117 gezelter 1490
118 gezelter 1930 //calculate atoms in cutoff groups
119     int nAtomsInGroups = 0;
120     int nCutoffGroupsInStamp = molStamp->getNCutoffGroups();
121    
122     for (int j=0; j < nCutoffGroupsInStamp; j++) {
123 tim 2469 cgStamp = molStamp->getCutoffGroupStamp(j);
124 gezelter 2204 nAtomsInGroups += cgStamp->getNMembers();
125 gezelter 1930 }
126 gezelter 1490
127 gezelter 1930 nGroups += nCutoffGroupsInStamp * nMolWithSameStamp;
128 chrisfen 2344
129 gezelter 1930 nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp;
130 gezelter 1490
131 gezelter 1930 //calculate atoms in rigid bodies
132     int nAtomsInRigidBodies = 0;
133 tim 1958 int nRigidBodiesInStamp = molStamp->getNRigidBodies();
134 gezelter 1930
135     for (int j=0; j < nRigidBodiesInStamp; j++) {
136 tim 2469 rbStamp = molStamp->getRigidBodyStamp(j);
137 gezelter 2204 nAtomsInRigidBodies += rbStamp->getNMembers();
138 gezelter 1930 }
139 gezelter 1490
140 gezelter 1930 nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp;
141     nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp;
142    
143 gezelter 2204 }
144 chrisfen 1636
145 chrisfen 2344 //every free atom (atom does not belong to cutoff groups) is a cutoff
146     //group therefore the total number of cutoff groups in the system is
147     //equal to the total number of atoms minus number of atoms belong to
148     //cutoff group defined in meta-data file plus the number of cutoff
149     //groups defined in meta-data file
150 gezelter 2204 nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups;
151 gezelter 1490
152 chrisfen 2344 //every free atom (atom does not belong to rigid bodies) is an
153     //integrable object therefore the total number of integrable objects
154     //in the system is equal to the total number of atoms minus number of
155     //atoms belong to rigid body defined in meta-data file plus the number
156     //of rigid bodies defined in meta-data file
157     nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms
158     + nGlobalRigidBodies_;
159    
160 gezelter 2204 nGlobalMols_ = molStampIds_.size();
161     molToProcMap_.resize(nGlobalMols_);
162     }
163 gezelter 1490
164 gezelter 2204 SimInfo::~SimInfo() {
165 tim 2082 std::map<int, Molecule*>::iterator i;
166     for (i = molecules_.begin(); i != molecules_.end(); ++i) {
167 gezelter 2204 delete i->second;
168 tim 2082 }
169     molecules_.clear();
170 tim 2187
171 gezelter 1930 delete sman_;
172     delete simParams_;
173     delete forceField_;
174 gezelter 2204 }
175 gezelter 1490
176 gezelter 2204 int SimInfo::getNGlobalConstraints() {
177 gezelter 1930 int nGlobalConstraints;
178     #ifdef IS_MPI
179     MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM,
180     MPI_COMM_WORLD);
181     #else
182     nGlobalConstraints = nConstraints_;
183     #endif
184     return nGlobalConstraints;
185 gezelter 2204 }
186 gezelter 1490
187 gezelter 2204 bool SimInfo::addMolecule(Molecule* mol) {
188 gezelter 1930 MoleculeIterator i;
189 gezelter 1490
190 gezelter 1930 i = molecules_.find(mol->getGlobalIndex());
191     if (i == molecules_.end() ) {
192 gezelter 1490
193 gezelter 2204 molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol));
194 gezelter 1930
195 gezelter 2204 nAtoms_ += mol->getNAtoms();
196     nBonds_ += mol->getNBonds();
197     nBends_ += mol->getNBends();
198     nTorsions_ += mol->getNTorsions();
199 gezelter 3432 nInversions_ += mol->getNInversions();
200 gezelter 2204 nRigidBodies_ += mol->getNRigidBodies();
201     nIntegrableObjects_ += mol->getNIntegrableObjects();
202     nCutoffGroups_ += mol->getNCutoffGroups();
203     nConstraints_ += mol->getNConstraintPairs();
204 gezelter 1490
205 gezelter 3442 addInteractionPairs(mol);
206    
207 gezelter 2204 return true;
208 gezelter 1930 } else {
209 gezelter 2204 return false;
210 gezelter 1930 }
211 gezelter 2204 }
212 gezelter 1490
213 gezelter 2204 bool SimInfo::removeMolecule(Molecule* mol) {
214 gezelter 1930 MoleculeIterator i;
215     i = molecules_.find(mol->getGlobalIndex());
216 gezelter 1490
217 gezelter 1930 if (i != molecules_.end() ) {
218 gezelter 1490
219 gezelter 2204 assert(mol == i->second);
220 gezelter 1930
221 gezelter 2204 nAtoms_ -= mol->getNAtoms();
222     nBonds_ -= mol->getNBonds();
223     nBends_ -= mol->getNBends();
224     nTorsions_ -= mol->getNTorsions();
225 gezelter 3432 nInversions_ -= mol->getNInversions();
226 gezelter 2204 nRigidBodies_ -= mol->getNRigidBodies();
227     nIntegrableObjects_ -= mol->getNIntegrableObjects();
228     nCutoffGroups_ -= mol->getNCutoffGroups();
229     nConstraints_ -= mol->getNConstraintPairs();
230 gezelter 1490
231 gezelter 3442 removeInteractionPairs(mol);
232 gezelter 2204 molecules_.erase(mol->getGlobalIndex());
233 gezelter 1490
234 gezelter 2204 delete mol;
235 gezelter 1930
236 gezelter 2204 return true;
237 gezelter 1930 } else {
238 gezelter 2204 return false;
239 gezelter 1930 }
240    
241    
242 gezelter 2204 }
243 gezelter 1930
244    
245 gezelter 2204 Molecule* SimInfo::beginMolecule(MoleculeIterator& i) {
246 gezelter 1930 i = molecules_.begin();
247     return i == molecules_.end() ? NULL : i->second;
248 gezelter 2204 }
249 gezelter 1930
250 gezelter 2204 Molecule* SimInfo::nextMolecule(MoleculeIterator& i) {
251 gezelter 1930 ++i;
252     return i == molecules_.end() ? NULL : i->second;
253 gezelter 2204 }
254 gezelter 1490
255    
256 gezelter 2204 void SimInfo::calcNdf() {
257 gezelter 1930 int ndf_local;
258     MoleculeIterator i;
259     std::vector<StuntDouble*>::iterator j;
260     Molecule* mol;
261     StuntDouble* integrableObject;
262 gezelter 1490
263 gezelter 1930 ndf_local = 0;
264    
265     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
266 gezelter 2204 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
267     integrableObject = mol->nextIntegrableObject(j)) {
268 gezelter 1490
269 gezelter 2204 ndf_local += 3;
270 gezelter 1490
271 gezelter 2204 if (integrableObject->isDirectional()) {
272     if (integrableObject->isLinear()) {
273     ndf_local += 2;
274     } else {
275     ndf_local += 3;
276     }
277     }
278 gezelter 1930
279 tim 2469 }
280     }
281 gezelter 1930
282     // n_constraints is local, so subtract them on each processor
283     ndf_local -= nConstraints_;
284    
285     #ifdef IS_MPI
286     MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
287     #else
288     ndf_ = ndf_local;
289     #endif
290    
291     // nZconstraints_ is global, as are the 3 COM translations for the
292     // entire system:
293     ndf_ = ndf_ - 3 - nZconstraint_;
294    
295 gezelter 2204 }
296 gezelter 1490
297 gezelter 2733 int SimInfo::getFdf() {
298     #ifdef IS_MPI
299     MPI_Allreduce(&fdf_local,&fdf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
300     #else
301     fdf_ = fdf_local;
302     #endif
303     return fdf_;
304     }
305    
306 gezelter 2204 void SimInfo::calcNdfRaw() {
307 gezelter 1930 int ndfRaw_local;
308 gezelter 1490
309 gezelter 1930 MoleculeIterator i;
310     std::vector<StuntDouble*>::iterator j;
311     Molecule* mol;
312     StuntDouble* integrableObject;
313    
314     // Raw degrees of freedom that we have to set
315     ndfRaw_local = 0;
316    
317     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
318 gezelter 2204 for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
319     integrableObject = mol->nextIntegrableObject(j)) {
320 gezelter 1930
321 gezelter 2204 ndfRaw_local += 3;
322 gezelter 1930
323 gezelter 2204 if (integrableObject->isDirectional()) {
324     if (integrableObject->isLinear()) {
325     ndfRaw_local += 2;
326     } else {
327     ndfRaw_local += 3;
328     }
329     }
330 gezelter 1930
331 gezelter 2204 }
332 gezelter 1930 }
333    
334     #ifdef IS_MPI
335     MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
336     #else
337     ndfRaw_ = ndfRaw_local;
338     #endif
339 gezelter 2204 }
340 gezelter 1490
341 gezelter 2204 void SimInfo::calcNdfTrans() {
342 gezelter 1930 int ndfTrans_local;
343 gezelter 1490
344 gezelter 1930 ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_;
345 gezelter 1490
346    
347 gezelter 1930 #ifdef IS_MPI
348     MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD);
349     #else
350     ndfTrans_ = ndfTrans_local;
351     #endif
352 gezelter 1490
353 gezelter 1930 ndfTrans_ = ndfTrans_ - 3 - nZconstraint_;
354    
355 gezelter 2204 }
356 gezelter 1490
357 gezelter 3442 void SimInfo::addInteractionPairs(Molecule* mol) {
358     ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
359 gezelter 1930 std::vector<Bond*>::iterator bondIter;
360     std::vector<Bend*>::iterator bendIter;
361     std::vector<Torsion*>::iterator torsionIter;
362 gezelter 3432 std::vector<Inversion*>::iterator inversionIter;
363 gezelter 1930 Bond* bond;
364     Bend* bend;
365     Torsion* torsion;
366 gezelter 3432 Inversion* inversion;
367 gezelter 1930 int a;
368     int b;
369     int c;
370     int d;
371 tim 2448
372 gezelter 3442 // atomGroups can be used to add special interaction maps between
373     // groups of atoms that are in two separate rigid bodies.
374     // However, most site-site interactions between two rigid bodies
375     // are probably not special, just the ones between the physically
376     // bonded atoms. Interactions *within* a single rigid body should
377     // always be excluded. These are done at the bottom of this
378     // function.
379    
380 tim 2448 std::map<int, std::set<int> > atomGroups;
381     Molecule::RigidBodyIterator rbIter;
382     RigidBody* rb;
383     Molecule::IntegrableObjectIterator ii;
384     StuntDouble* integrableObject;
385 gezelter 1930
386 gezelter 3442 for (integrableObject = mol->beginIntegrableObject(ii);
387     integrableObject != NULL;
388     integrableObject = mol->nextIntegrableObject(ii)) {
389    
390 tim 2448 if (integrableObject->isRigidBody()) {
391 gezelter 3442 rb = static_cast<RigidBody*>(integrableObject);
392     std::vector<Atom*> atoms = rb->getAtoms();
393     std::set<int> rigidAtoms;
394     for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
395     rigidAtoms.insert(atoms[i]->getGlobalIndex());
396     }
397     for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
398     atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
399     }
400 tim 2448 } else {
401     std::set<int> oneAtomSet;
402     oneAtomSet.insert(integrableObject->getGlobalIndex());
403     atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));
404     }
405     }
406 gezelter 3442
407     for (bond= mol->beginBond(bondIter); bond != NULL;
408     bond = mol->nextBond(bondIter)) {
409 tim 2448
410 gezelter 3442 a = bond->getAtomA()->getGlobalIndex();
411     b = bond->getAtomB()->getGlobalIndex();
412 tim 2448
413 gezelter 3442 if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
414     oneTwoInteractions_.addPair(a, b);
415     } else {
416     excludedInteractions_.addPair(a, b);
417     }
418 gezelter 1930 }
419 gezelter 1490
420 gezelter 3442 for (bend= mol->beginBend(bendIter); bend != NULL;
421     bend = mol->nextBend(bendIter)) {
422    
423 gezelter 2204 a = bend->getAtomA()->getGlobalIndex();
424     b = bend->getAtomB()->getGlobalIndex();
425     c = bend->getAtomC()->getGlobalIndex();
426 gezelter 3442
427     if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
428     oneTwoInteractions_.addPair(a, b);
429     oneTwoInteractions_.addPair(b, c);
430     } else {
431     excludedInteractions_.addPair(a, b);
432     excludedInteractions_.addPair(b, c);
433     }
434 gezelter 1490
435 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
436     oneThreeInteractions_.addPair(a, c);
437     } else {
438     excludedInteractions_.addPair(a, c);
439     }
440 gezelter 1930 }
441 gezelter 1490
442 gezelter 3442 for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
443     torsion = mol->nextTorsion(torsionIter)) {
444    
445 gezelter 2204 a = torsion->getAtomA()->getGlobalIndex();
446     b = torsion->getAtomB()->getGlobalIndex();
447     c = torsion->getAtomC()->getGlobalIndex();
448 gezelter 3442 d = torsion->getAtomD()->getGlobalIndex();
449 cli2 3446
450 gezelter 3442 if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
451     oneTwoInteractions_.addPair(a, b);
452     oneTwoInteractions_.addPair(b, c);
453     oneTwoInteractions_.addPair(c, d);
454     } else {
455     excludedInteractions_.addPair(a, b);
456     excludedInteractions_.addPair(b, c);
457     excludedInteractions_.addPair(c, d);
458     }
459 gezelter 1490
460 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
461     oneThreeInteractions_.addPair(a, c);
462     oneThreeInteractions_.addPair(b, d);
463     } else {
464     excludedInteractions_.addPair(a, c);
465     excludedInteractions_.addPair(b, d);
466     }
467 tim 2448
468 gezelter 3442 if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
469     oneFourInteractions_.addPair(a, d);
470     } else {
471     excludedInteractions_.addPair(a, d);
472     }
473 gezelter 1490 }
474    
475 gezelter 3432 for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
476     inversion = mol->nextInversion(inversionIter)) {
477 gezelter 3442
478 gezelter 3432 a = inversion->getAtomA()->getGlobalIndex();
479     b = inversion->getAtomB()->getGlobalIndex();
480     c = inversion->getAtomC()->getGlobalIndex();
481     d = inversion->getAtomD()->getGlobalIndex();
482    
483 gezelter 3442 if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
484     oneTwoInteractions_.addPair(a, b);
485     oneTwoInteractions_.addPair(a, c);
486     oneTwoInteractions_.addPair(a, d);
487     } else {
488     excludedInteractions_.addPair(a, b);
489     excludedInteractions_.addPair(a, c);
490     excludedInteractions_.addPair(a, d);
491     }
492 gezelter 3432
493 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
494     oneThreeInteractions_.addPair(b, c);
495     oneThreeInteractions_.addPair(b, d);
496     oneThreeInteractions_.addPair(c, d);
497     } else {
498     excludedInteractions_.addPair(b, c);
499     excludedInteractions_.addPair(b, d);
500     excludedInteractions_.addPair(c, d);
501     }
502 gezelter 3432 }
503    
504 gezelter 3442 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
505     rb = mol->nextRigidBody(rbIter)) {
506 gezelter 2204 std::vector<Atom*> atoms = rb->getAtoms();
507 gezelter 3442 for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
508     for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
509 gezelter 2204 a = atoms[i]->getGlobalIndex();
510     b = atoms[j]->getGlobalIndex();
511 gezelter 3442 excludedInteractions_.addPair(a, b);
512 gezelter 2204 }
513     }
514 tim 2114 }
515    
516 gezelter 2204 }
517 gezelter 1930
518 gezelter 3442 void SimInfo::removeInteractionPairs(Molecule* mol) {
519     ForceFieldOptions& options_ = forceField_->getForceFieldOptions();
520 gezelter 1930 std::vector<Bond*>::iterator bondIter;
521     std::vector<Bend*>::iterator bendIter;
522     std::vector<Torsion*>::iterator torsionIter;
523 gezelter 3432 std::vector<Inversion*>::iterator inversionIter;
524 gezelter 1930 Bond* bond;
525     Bend* bend;
526     Torsion* torsion;
527 gezelter 3432 Inversion* inversion;
528 gezelter 1930 int a;
529     int b;
530     int c;
531     int d;
532 tim 2448
533     std::map<int, std::set<int> > atomGroups;
534     Molecule::RigidBodyIterator rbIter;
535     RigidBody* rb;
536     Molecule::IntegrableObjectIterator ii;
537     StuntDouble* integrableObject;
538 gezelter 1930
539 gezelter 3442 for (integrableObject = mol->beginIntegrableObject(ii);
540     integrableObject != NULL;
541     integrableObject = mol->nextIntegrableObject(ii)) {
542    
543 tim 2448 if (integrableObject->isRigidBody()) {
544 gezelter 3442 rb = static_cast<RigidBody*>(integrableObject);
545     std::vector<Atom*> atoms = rb->getAtoms();
546     std::set<int> rigidAtoms;
547     for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
548     rigidAtoms.insert(atoms[i]->getGlobalIndex());
549     }
550     for (int i = 0; i < static_cast<int>(atoms.size()); ++i) {
551     atomGroups.insert(std::map<int, std::set<int> >::value_type(atoms[i]->getGlobalIndex(), rigidAtoms));
552     }
553 tim 2448 } else {
554     std::set<int> oneAtomSet;
555     oneAtomSet.insert(integrableObject->getGlobalIndex());
556     atomGroups.insert(std::map<int, std::set<int> >::value_type(integrableObject->getGlobalIndex(), oneAtomSet));
557     }
558     }
559    
560 gezelter 3442 for (bond= mol->beginBond(bondIter); bond != NULL;
561     bond = mol->nextBond(bondIter)) {
562    
563     a = bond->getAtomA()->getGlobalIndex();
564     b = bond->getAtomB()->getGlobalIndex();
565 tim 2448
566 gezelter 3442 if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
567     oneTwoInteractions_.removePair(a, b);
568     } else {
569     excludedInteractions_.removePair(a, b);
570     }
571 gezelter 1490 }
572 gezelter 1930
573 gezelter 3442 for (bend= mol->beginBend(bendIter); bend != NULL;
574     bend = mol->nextBend(bendIter)) {
575    
576 gezelter 2204 a = bend->getAtomA()->getGlobalIndex();
577     b = bend->getAtomB()->getGlobalIndex();
578     c = bend->getAtomC()->getGlobalIndex();
579 gezelter 3442
580     if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
581     oneTwoInteractions_.removePair(a, b);
582     oneTwoInteractions_.removePair(b, c);
583     } else {
584     excludedInteractions_.removePair(a, b);
585     excludedInteractions_.removePair(b, c);
586     }
587 gezelter 1930
588 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
589     oneThreeInteractions_.removePair(a, c);
590     } else {
591     excludedInteractions_.removePair(a, c);
592     }
593 gezelter 1490 }
594 gezelter 1930
595 gezelter 3442 for (torsion= mol->beginTorsion(torsionIter); torsion != NULL;
596     torsion = mol->nextTorsion(torsionIter)) {
597    
598 gezelter 2204 a = torsion->getAtomA()->getGlobalIndex();
599     b = torsion->getAtomB()->getGlobalIndex();
600     c = torsion->getAtomC()->getGlobalIndex();
601 gezelter 3442 d = torsion->getAtomD()->getGlobalIndex();
602    
603     if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
604     oneTwoInteractions_.removePair(a, b);
605     oneTwoInteractions_.removePair(b, c);
606     oneTwoInteractions_.removePair(c, d);
607     } else {
608     excludedInteractions_.removePair(a, b);
609     excludedInteractions_.removePair(b, c);
610     excludedInteractions_.removePair(c, d);
611     }
612 gezelter 1930
613 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
614     oneThreeInteractions_.removePair(a, c);
615     oneThreeInteractions_.removePair(b, d);
616     } else {
617     excludedInteractions_.removePair(a, c);
618     excludedInteractions_.removePair(b, d);
619     }
620 tim 2448
621 gezelter 3442 if (options_.havevdw14scale() || options_.haveelectrostatic14scale()) {
622     oneFourInteractions_.removePair(a, d);
623     } else {
624     excludedInteractions_.removePair(a, d);
625     }
626     }
627 tim 2448
628 gezelter 3442 for (inversion= mol->beginInversion(inversionIter); inversion != NULL;
629     inversion = mol->nextInversion(inversionIter)) {
630 tim 2448
631 gezelter 3432 a = inversion->getAtomA()->getGlobalIndex();
632     b = inversion->getAtomB()->getGlobalIndex();
633     c = inversion->getAtomC()->getGlobalIndex();
634     d = inversion->getAtomD()->getGlobalIndex();
635    
636 gezelter 3442 if (options_.havevdw12scale() || options_.haveelectrostatic12scale()) {
637     oneTwoInteractions_.removePair(a, b);
638     oneTwoInteractions_.removePair(a, c);
639     oneTwoInteractions_.removePair(a, d);
640     } else {
641     excludedInteractions_.removePair(a, b);
642     excludedInteractions_.removePair(a, c);
643     excludedInteractions_.removePair(a, d);
644     }
645 gezelter 3432
646 gezelter 3442 if (options_.havevdw13scale() || options_.haveelectrostatic13scale()) {
647     oneThreeInteractions_.removePair(b, c);
648     oneThreeInteractions_.removePair(b, d);
649     oneThreeInteractions_.removePair(c, d);
650     } else {
651     excludedInteractions_.removePair(b, c);
652     excludedInteractions_.removePair(b, d);
653     excludedInteractions_.removePair(c, d);
654     }
655 gezelter 3432 }
656    
657 gezelter 3442 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
658     rb = mol->nextRigidBody(rbIter)) {
659 gezelter 2204 std::vector<Atom*> atoms = rb->getAtoms();
660 gezelter 3442 for (int i = 0; i < static_cast<int>(atoms.size()) -1 ; ++i) {
661     for (int j = i + 1; j < static_cast<int>(atoms.size()); ++j) {
662 gezelter 2204 a = atoms[i]->getGlobalIndex();
663     b = atoms[j]->getGlobalIndex();
664 gezelter 3442 excludedInteractions_.removePair(a, b);
665 gezelter 2204 }
666     }
667 tim 2114 }
668 gezelter 3442
669 gezelter 2204 }
670 gezelter 3442
671    
672 gezelter 2204 void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) {
673 gezelter 1930 int curStampId;
674 gezelter 3442
675 gezelter 1930 //index from 0
676     curStampId = moleculeStamps_.size();
677 gezelter 1490
678 gezelter 1930 moleculeStamps_.push_back(molStamp);
679     molStampIds_.insert(molStampIds_.end(), nmol, curStampId);
680 gezelter 2204 }
681 gezelter 1490
682 gezelter 2204 void SimInfo::update() {
683 gezelter 1490
684 gezelter 1930 setupSimType();
685 gezelter 1490
686 gezelter 1930 #ifdef IS_MPI
687     setupFortranParallel();
688     #endif
689 gezelter 1490
690 gezelter 1930 setupFortranSim();
691 gezelter 1490
692 gezelter 1930 //setup fortran force field
693     /** @deprecate */
694     int isError = 0;
695 chrisfen 2297
696 chrisfen 3013 setupCutoff();
697    
698 chrisfen 2302 setupElectrostaticSummationMethod( isError );
699 chrisfen 2425 setupSwitchingFunction();
700 chrisfen 2917 setupAccumulateBoxDipole();
701 chrisfen 2297
702 gezelter 1930 if(isError){
703 gezelter 2204 sprintf( painCave.errMsg,
704     "ForceField error: There was an error initializing the forceField in fortran.\n" );
705     painCave.isFatal = 1;
706     simError();
707 gezelter 1930 }
708 gezelter 1490
709 gezelter 1930 calcNdf();
710     calcNdfRaw();
711     calcNdfTrans();
712    
713     fortranInitialized_ = true;
714 gezelter 2204 }
715 gezelter 1490
716 gezelter 2204 std::set<AtomType*> SimInfo::getUniqueAtomTypes() {
717 gezelter 1930 SimInfo::MoleculeIterator mi;
718     Molecule* mol;
719     Molecule::AtomIterator ai;
720     Atom* atom;
721     std::set<AtomType*> atomTypes;
722 gezelter 1490
723 gezelter 1930 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
724 gezelter 1490
725 gezelter 2204 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
726     atomTypes.insert(atom->getAtomType());
727     }
728 gezelter 1930
729     }
730 gezelter 1490
731 gezelter 1930 return atomTypes;
732 gezelter 2204 }
733 gezelter 1490
734 gezelter 2204 void SimInfo::setupSimType() {
735 gezelter 1930 std::set<AtomType*>::iterator i;
736     std::set<AtomType*> atomTypes;
737     atomTypes = getUniqueAtomTypes();
738 gezelter 1490
739 gezelter 1930 int useLennardJones = 0;
740     int useElectrostatic = 0;
741     int useEAM = 0;
742 chuckv 2433 int useSC = 0;
743 gezelter 1930 int useCharge = 0;
744     int useDirectional = 0;
745     int useDipole = 0;
746     int useGayBerne = 0;
747     int useSticky = 0;
748 chrisfen 2220 int useStickyPower = 0;
749 gezelter 1930 int useShape = 0;
750     int useFLARB = 0; //it is not in AtomType yet
751     int useDirectionalAtom = 0;
752     int useElectrostatics = 0;
753     //usePBC and useRF are from simParams
754 tim 2364 int usePBC = simParams_->getUsePeriodicBoundaryConditions();
755 chrisfen 2310 int useRF;
756 chrisfen 2419 int useSF;
757 chrisfen 2917 int useSP;
758     int useBoxDipole;
759 gezelter 3126
760 tim 2364 std::string myMethod;
761 gezelter 1490
762 chrisfen 2310 // set the useRF logical
763 tim 2364 useRF = 0;
764 chrisfen 2419 useSF = 0;
765 gezelter 3054 useSP = 0;
766 gezelter 3470 useBoxDipole = 0;
767 chrisfen 2390
768    
769 tim 2364 if (simParams_->haveElectrostaticSummationMethod()) {
770 chrisfen 2390 std::string myMethod = simParams_->getElectrostaticSummationMethod();
771     toUpper(myMethod);
772 chrisfen 2917 if (myMethod == "REACTION_FIELD"){
773 gezelter 3054 useRF = 1;
774 chrisfen 2917 } else if (myMethod == "SHIFTED_FORCE"){
775     useSF = 1;
776     } else if (myMethod == "SHIFTED_POTENTIAL"){
777     useSP = 1;
778 chrisfen 2390 }
779 tim 2364 }
780 chrisfen 2917
781     if (simParams_->haveAccumulateBoxDipole())
782     if (simParams_->getAccumulateBoxDipole())
783     useBoxDipole = 1;
784 chrisfen 2310
785 gezelter 3126 useAtomicVirial_ = simParams_->getUseAtomicVirial();
786    
787 gezelter 1930 //loop over all of the atom types
788     for (i = atomTypes.begin(); i != atomTypes.end(); ++i) {
789 gezelter 2204 useLennardJones |= (*i)->isLennardJones();
790     useElectrostatic |= (*i)->isElectrostatic();
791     useEAM |= (*i)->isEAM();
792 chuckv 2433 useSC |= (*i)->isSC();
793 gezelter 2204 useCharge |= (*i)->isCharge();
794     useDirectional |= (*i)->isDirectional();
795     useDipole |= (*i)->isDipole();
796     useGayBerne |= (*i)->isGayBerne();
797     useSticky |= (*i)->isSticky();
798 chrisfen 2220 useStickyPower |= (*i)->isStickyPower();
799 gezelter 2204 useShape |= (*i)->isShape();
800 gezelter 1930 }
801 gezelter 1490
802 chrisfen 2220 if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) {
803 gezelter 2204 useDirectionalAtom = 1;
804 gezelter 1930 }
805 gezelter 1490
806 gezelter 1930 if (useCharge || useDipole) {
807 gezelter 2204 useElectrostatics = 1;
808 gezelter 1930 }
809 gezelter 1490
810 gezelter 1930 #ifdef IS_MPI
811     int temp;
812 gezelter 1490
813 gezelter 1930 temp = usePBC;
814     MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
815 gezelter 1490
816 gezelter 1930 temp = useDirectionalAtom;
817     MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
818 gezelter 1490
819 gezelter 1930 temp = useLennardJones;
820     MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
821 gezelter 1490
822 gezelter 1930 temp = useElectrostatics;
823     MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
824 gezelter 1490
825 gezelter 1930 temp = useCharge;
826     MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
827 gezelter 1490
828 gezelter 1930 temp = useDipole;
829     MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
830 gezelter 1490
831 gezelter 1930 temp = useSticky;
832     MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
833 gezelter 1490
834 chrisfen 2220 temp = useStickyPower;
835     MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
836    
837 gezelter 1930 temp = useGayBerne;
838     MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
839 gezelter 1490
840 gezelter 1930 temp = useEAM;
841     MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
842 gezelter 1490
843 chuckv 2433 temp = useSC;
844     MPI_Allreduce(&temp, &useSC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
845    
846 gezelter 1930 temp = useShape;
847     MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
848    
849     temp = useFLARB;
850     MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
851    
852 chrisfen 2310 temp = useRF;
853     MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
854    
855 chrisfen 2419 temp = useSF;
856 chrisfen 2917 MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
857 chrisfen 2404
858 chrisfen 2917 temp = useSP;
859     MPI_Allreduce(&temp, &useSP, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
860    
861     temp = useBoxDipole;
862     MPI_Allreduce(&temp, &useBoxDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
863    
864 gezelter 3126 temp = useAtomicVirial_;
865     MPI_Allreduce(&temp, &useAtomicVirial_, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD);
866    
867 gezelter 1490 #endif
868    
869 gezelter 1930 fInfo_.SIM_uses_PBC = usePBC;
870     fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom;
871     fInfo_.SIM_uses_LennardJones = useLennardJones;
872     fInfo_.SIM_uses_Electrostatics = useElectrostatics;
873     fInfo_.SIM_uses_Charges = useCharge;
874     fInfo_.SIM_uses_Dipoles = useDipole;
875     fInfo_.SIM_uses_Sticky = useSticky;
876 chrisfen 2220 fInfo_.SIM_uses_StickyPower = useStickyPower;
877 gezelter 1930 fInfo_.SIM_uses_GayBerne = useGayBerne;
878     fInfo_.SIM_uses_EAM = useEAM;
879 chuckv 2433 fInfo_.SIM_uses_SC = useSC;
880 gezelter 1930 fInfo_.SIM_uses_Shapes = useShape;
881     fInfo_.SIM_uses_FLARB = useFLARB;
882 chrisfen 2310 fInfo_.SIM_uses_RF = useRF;
883 chrisfen 2419 fInfo_.SIM_uses_SF = useSF;
884 chrisfen 2917 fInfo_.SIM_uses_SP = useSP;
885     fInfo_.SIM_uses_BoxDipole = useBoxDipole;
886 gezelter 3126 fInfo_.SIM_uses_AtomicVirial = useAtomicVirial_;
887 gezelter 2204 }
888 gezelter 1490
889 gezelter 2204 void SimInfo::setupFortranSim() {
890 gezelter 1930 int isError;
891 gezelter 3442 int nExclude, nOneTwo, nOneThree, nOneFour;
892 gezelter 1930 std::vector<int> fortranGlobalGroupMembership;
893    
894     isError = 0;
895 gezelter 1490
896 gezelter 1930 //globalGroupMembership_ is filled by SimCreator
897     for (int i = 0; i < nGlobalAtoms_; i++) {
898 gezelter 2204 fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1);
899 gezelter 1930 }
900 gezelter 1490
901 gezelter 1930 //calculate mass ratio of cutoff group
902 tim 2759 std::vector<RealType> mfact;
903 gezelter 1930 SimInfo::MoleculeIterator mi;
904     Molecule* mol;
905     Molecule::CutoffGroupIterator ci;
906     CutoffGroup* cg;
907     Molecule::AtomIterator ai;
908     Atom* atom;
909 tim 2759 RealType totalMass;
910 gezelter 1930
911     //to avoid memory reallocation, reserve enough space for mfact
912     mfact.reserve(getNCutoffGroups());
913 gezelter 1490
914 gezelter 1930 for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
915 gezelter 2204 for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
916 gezelter 1490
917 gezelter 2204 totalMass = cg->getMass();
918     for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) {
919 chrisfen 2344 // Check for massless groups - set mfact to 1 if true
920     if (totalMass != 0)
921     mfact.push_back(atom->getMass()/totalMass);
922     else
923     mfact.push_back( 1.0 );
924 gezelter 2204 }
925     }
926 gezelter 1930 }
927 gezelter 1490
928 gezelter 1930 //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!)
929     std::vector<int> identArray;
930 gezelter 1490
931 gezelter 1930 //to avoid memory reallocation, reserve enough space identArray
932     identArray.reserve(getNAtoms());
933    
934     for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
935 gezelter 2204 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
936     identArray.push_back(atom->getIdent());
937     }
938 gezelter 1930 }
939 gezelter 1490
940 gezelter 1930 //fill molMembershipArray
941     //molMembershipArray is filled by SimCreator
942     std::vector<int> molMembershipArray(nGlobalAtoms_);
943     for (int i = 0; i < nGlobalAtoms_; i++) {
944 gezelter 2204 molMembershipArray[i] = globalMolMembership_[i] + 1;
945 gezelter 1930 }
946    
947     //setup fortran simulation
948 gezelter 3442
949     nExclude = excludedInteractions_.getSize();
950     nOneTwo = oneTwoInteractions_.getSize();
951     nOneThree = oneThreeInteractions_.getSize();
952     nOneFour = oneFourInteractions_.getSize();
953    
954     int* excludeList = excludedInteractions_.getPairList();
955     int* oneTwoList = oneTwoInteractions_.getPairList();
956     int* oneThreeList = oneThreeInteractions_.getPairList();
957     int* oneFourList = oneFourInteractions_.getPairList();
958    
959 gezelter 3388 setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0],
960 gezelter 3442 &nExclude, excludeList,
961     &nOneTwo, oneTwoList,
962     &nOneThree, oneThreeList,
963     &nOneFour, oneFourList,
964 gezelter 3388 &molMembershipArray[0], &mfact[0], &nCutoffGroups_,
965     &fortranGlobalGroupMembership[0], &isError);
966    
967 gezelter 1930 if( isError ){
968 gezelter 3388
969 gezelter 2204 sprintf( painCave.errMsg,
970     "There was an error setting the simulation information in fortran.\n" );
971     painCave.isFatal = 1;
972     painCave.severity = OOPSE_ERROR;
973     simError();
974 gezelter 1930 }
975 gezelter 3388
976    
977 gezelter 1930 sprintf( checkPointMsg,
978 gezelter 2204 "succesfully sent the simulation information to fortran.\n");
979 gezelter 3388
980     errorCheckPoint();
981    
982 chuckv 3080 // Setup number of neighbors in neighbor list if present
983     if (simParams_->haveNeighborListNeighbors()) {
984 chuckv 3121 int nlistNeighbors = simParams_->getNeighborListNeighbors();
985     setNeighbors(&nlistNeighbors);
986 chuckv 3080 }
987    
988    
989 gezelter 2204 }
990 gezelter 1490
991    
992 gezelter 2204 void SimInfo::setupFortranParallel() {
993 gezelter 3388 #ifdef IS_MPI
994 gezelter 1930 //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex
995     std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0);
996     std::vector<int> localToGlobalCutoffGroupIndex;
997     SimInfo::MoleculeIterator mi;
998     Molecule::AtomIterator ai;
999     Molecule::CutoffGroupIterator ci;
1000     Molecule* mol;
1001     Atom* atom;
1002     CutoffGroup* cg;
1003     mpiSimData parallelData;
1004     int isError;
1005 gezelter 1490
1006 gezelter 1930 for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1007 gezelter 1490
1008 gezelter 2204 //local index(index in DataStorge) of atom is important
1009     for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
1010     localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1;
1011     }
1012 gezelter 1490
1013 gezelter 2204 //local index of cutoff group is trivial, it only depends on the order of travesing
1014     for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) {
1015     localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1);
1016     }
1017 gezelter 1930
1018     }
1019 gezelter 1490
1020 gezelter 1930 //fill up mpiSimData struct
1021     parallelData.nMolGlobal = getNGlobalMolecules();
1022     parallelData.nMolLocal = getNMolecules();
1023     parallelData.nAtomsGlobal = getNGlobalAtoms();
1024     parallelData.nAtomsLocal = getNAtoms();
1025     parallelData.nGroupsGlobal = getNGlobalCutoffGroups();
1026     parallelData.nGroupsLocal = getNCutoffGroups();
1027     parallelData.myNode = worldRank;
1028     MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors));
1029 gezelter 1490
1030 gezelter 1930 //pass mpiSimData struct and index arrays to fortran
1031     setFsimParallel(&parallelData, &(parallelData.nAtomsLocal),
1032     &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal),
1033     &localToGlobalCutoffGroupIndex[0], &isError);
1034 gezelter 1490
1035 gezelter 1930 if (isError) {
1036 gezelter 2204 sprintf(painCave.errMsg,
1037     "mpiRefresh errror: fortran didn't like something we gave it.\n");
1038     painCave.isFatal = 1;
1039     simError();
1040 gezelter 1930 }
1041 gezelter 1490
1042 gezelter 1930 sprintf(checkPointMsg, " mpiRefresh successful.\n");
1043 gezelter 3388 errorCheckPoint();
1044 gezelter 1490
1045 gezelter 3388 #endif
1046 gezelter 2204 }
1047 chrisfen 1636
1048 gezelter 2463 void SimInfo::setupCutoff() {
1049 gezelter 1490
1050 chuckv 2533 ForceFieldOptions& forceFieldOptions_ = forceField_->getForceFieldOptions();
1051    
1052 gezelter 2463 // Check the cutoff policy
1053 chuckv 2533 int cp = TRADITIONAL_CUTOFF_POLICY; // Set to traditional by default
1054    
1055 chrisfen 3129 // Set LJ shifting bools to false
1056 gezelter 3559 ljsp_ = 0;
1057     ljsf_ = 0;
1058 chrisfen 3129
1059 chuckv 2533 std::string myPolicy;
1060     if (forceFieldOptions_.haveCutoffPolicy()){
1061     myPolicy = forceFieldOptions_.getCutoffPolicy();
1062     }else if (simParams_->haveCutoffPolicy()) {
1063     myPolicy = simParams_->getCutoffPolicy();
1064     }
1065    
1066     if (!myPolicy.empty()){
1067 tim 2364 toUpper(myPolicy);
1068 gezelter 2285 if (myPolicy == "MIX") {
1069     cp = MIX_CUTOFF_POLICY;
1070     } else {
1071     if (myPolicy == "MAX") {
1072     cp = MAX_CUTOFF_POLICY;
1073     } else {
1074     if (myPolicy == "TRADITIONAL") {
1075     cp = TRADITIONAL_CUTOFF_POLICY;
1076     } else {
1077     // throw error
1078     sprintf( painCave.errMsg,
1079     "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() );
1080     painCave.isFatal = 1;
1081     simError();
1082     }
1083     }
1084     }
1085 gezelter 2463 }
1086     notifyFortranCutoffPolicy(&cp);
1087 chuckv 2328
1088 gezelter 2463 // Check the Skin Thickness for neighborlists
1089 tim 2759 RealType skin;
1090 gezelter 2463 if (simParams_->haveSkinThickness()) {
1091     skin = simParams_->getSkinThickness();
1092     notifyFortranSkinThickness(&skin);
1093     }
1094    
1095     // Check if the cutoff was set explicitly:
1096     if (simParams_->haveCutoffRadius()) {
1097     rcut_ = simParams_->getCutoffRadius();
1098     if (simParams_->haveSwitchingRadius()) {
1099     rsw_ = simParams_->getSwitchingRadius();
1100     } else {
1101 chrisfen 2578 if (fInfo_.SIM_uses_Charges |
1102     fInfo_.SIM_uses_Dipoles |
1103     fInfo_.SIM_uses_RF) {
1104    
1105     rsw_ = 0.85 * rcut_;
1106     sprintf(painCave.errMsg,
1107     "SimCreator Warning: No value was set for the switchingRadius.\n"
1108 chrisfen 2579 "\tOOPSE will use a default value of 85 percent of the cutoffRadius.\n"
1109 chrisfen 2578 "\tswitchingRadius = %f. for this simulation\n", rsw_);
1110     painCave.isFatal = 0;
1111     simError();
1112     } else {
1113     rsw_ = rcut_;
1114     sprintf(painCave.errMsg,
1115     "SimCreator Warning: No value was set for the switchingRadius.\n"
1116     "\tOOPSE will use the same value as the cutoffRadius.\n"
1117     "\tswitchingRadius = %f. for this simulation\n", rsw_);
1118     painCave.isFatal = 0;
1119     simError();
1120     }
1121 chrisfen 2579 }
1122 chrisfen 3129
1123     if (simParams_->haveElectrostaticSummationMethod()) {
1124     std::string myMethod = simParams_->getElectrostaticSummationMethod();
1125     toUpper(myMethod);
1126    
1127     if (myMethod == "SHIFTED_POTENTIAL") {
1128 gezelter 3559 ljsp_ = 1;
1129 chrisfen 3129 } else if (myMethod == "SHIFTED_FORCE") {
1130 gezelter 3559 ljsf_ = 1;
1131 chrisfen 3129 }
1132     }
1133 gezelter 3559
1134 chrisfen 3129 notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1135 chrisfen 2579
1136 gezelter 2463 } else {
1137    
1138     // For electrostatic atoms, we'll assume a large safe value:
1139     if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) {
1140     sprintf(painCave.errMsg,
1141     "SimCreator Warning: No value was set for the cutoffRadius.\n"
1142     "\tOOPSE will use a default value of 15.0 angstroms"
1143     "\tfor the cutoffRadius.\n");
1144     painCave.isFatal = 0;
1145     simError();
1146     rcut_ = 15.0;
1147    
1148     if (simParams_->haveElectrostaticSummationMethod()) {
1149     std::string myMethod = simParams_->getElectrostaticSummationMethod();
1150     toUpper(myMethod);
1151 chrisfen 3129
1152     // For the time being, we're tethering the LJ shifted behavior to the
1153     // electrostaticSummationMethod keyword options
1154     if (myMethod == "SHIFTED_POTENTIAL") {
1155 gezelter 3559 ljsp_ = 1;
1156 chrisfen 3129 } else if (myMethod == "SHIFTED_FORCE") {
1157 gezelter 3559 ljsf_ = 1;
1158 chrisfen 3129 }
1159     if (myMethod == "SHIFTED_POTENTIAL" || myMethod == "SHIFTED_FORCE") {
1160 gezelter 2463 if (simParams_->haveSwitchingRadius()){
1161     sprintf(painCave.errMsg,
1162     "SimInfo Warning: A value was set for the switchingRadius\n"
1163     "\teven though the electrostaticSummationMethod was\n"
1164     "\tset to %s\n", myMethod.c_str());
1165     painCave.isFatal = 1;
1166     simError();
1167     }
1168     }
1169     }
1170    
1171     if (simParams_->haveSwitchingRadius()){
1172     rsw_ = simParams_->getSwitchingRadius();
1173     } else {
1174     sprintf(painCave.errMsg,
1175     "SimCreator Warning: No value was set for switchingRadius.\n"
1176     "\tOOPSE will use a default value of\n"
1177     "\t0.85 * cutoffRadius for the switchingRadius\n");
1178     painCave.isFatal = 0;
1179     simError();
1180     rsw_ = 0.85 * rcut_;
1181     }
1182 chrisfen 3129
1183     notifyFortranCutoffs(&rcut_, &rsw_, &ljsp_, &ljsf_);
1184    
1185 gezelter 2463 } else {
1186     // We didn't set rcut explicitly, and we don't have electrostatic atoms, so
1187     // We'll punt and let fortran figure out the cutoffs later.
1188    
1189     notifyFortranYouAreOnYourOwn();
1190 chuckv 2328
1191 gezelter 2463 }
1192 chuckv 2328 }
1193 gezelter 2204 }
1194 gezelter 1490
1195 chrisfen 2302 void SimInfo::setupElectrostaticSummationMethod( int isError ) {
1196 chrisfen 2297
1197     int errorOut;
1198 chrisfen 2302 int esm = NONE;
1199 chrisfen 2408 int sm = UNDAMPED;
1200 tim 2759 RealType alphaVal;
1201     RealType dielectric;
1202 chrisfen 3013
1203 chrisfen 2297 errorOut = isError;
1204    
1205 chrisfen 2302 if (simParams_->haveElectrostaticSummationMethod()) {
1206 chrisfen 2303 std::string myMethod = simParams_->getElectrostaticSummationMethod();
1207 tim 2364 toUpper(myMethod);
1208 chrisfen 2302 if (myMethod == "NONE") {
1209     esm = NONE;
1210 chrisfen 2297 } else {
1211 chrisfen 2408 if (myMethod == "SWITCHING_FUNCTION") {
1212     esm = SWITCHING_FUNCTION;
1213 chrisfen 2297 } else {
1214 chrisfen 2408 if (myMethod == "SHIFTED_POTENTIAL") {
1215     esm = SHIFTED_POTENTIAL;
1216     } else {
1217     if (myMethod == "SHIFTED_FORCE") {
1218     esm = SHIFTED_FORCE;
1219 chrisfen 2297 } else {
1220 chrisfen 3020 if (myMethod == "REACTION_FIELD") {
1221 chrisfen 2408 esm = REACTION_FIELD;
1222 chrisfen 3020 dielectric = simParams_->getDielectric();
1223     if (!simParams_->haveDielectric()) {
1224     // throw warning
1225     sprintf( painCave.errMsg,
1226     "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
1227     "\tA default value of %f will be used for the dielectric.\n", dielectric);
1228     painCave.isFatal = 0;
1229     simError();
1230     }
1231 chrisfen 2408 } else {
1232     // throw error
1233     sprintf( painCave.errMsg,
1234 gezelter 2463 "SimInfo error: Unknown electrostaticSummationMethod.\n"
1235     "\t(Input file specified %s .)\n"
1236     "\telectrostaticSummationMethod must be one of: \"none\",\n"
1237     "\t\"shifted_potential\", \"shifted_force\", or \n"
1238     "\t\"reaction_field\".\n", myMethod.c_str() );
1239 chrisfen 2408 painCave.isFatal = 1;
1240     simError();
1241     }
1242     }
1243     }
1244 chrisfen 2297 }
1245     }
1246     }
1247 chrisfen 2408
1248 chrisfen 2415 if (simParams_->haveElectrostaticScreeningMethod()) {
1249     std::string myScreen = simParams_->getElectrostaticScreeningMethod();
1250 chrisfen 2408 toUpper(myScreen);
1251     if (myScreen == "UNDAMPED") {
1252     sm = UNDAMPED;
1253     } else {
1254     if (myScreen == "DAMPED") {
1255     sm = DAMPED;
1256     if (!simParams_->haveDampingAlpha()) {
1257 chrisfen 3013 // first set a cutoff dependent alpha value
1258     // we assume alpha depends linearly with rcut from 0 to 20.5 ang
1259     alphaVal = 0.5125 - rcut_* 0.025;
1260     // for values rcut > 20.5, alpha is zero
1261     if (alphaVal < 0) alphaVal = 0;
1262    
1263     // throw warning
1264 chrisfen 2408 sprintf( painCave.errMsg,
1265 gezelter 2463 "SimInfo warning: dampingAlpha was not specified in the input file.\n"
1266 chrisfen 3013 "\tA default value of %f (1/ang) will be used for the cutoff of\n\t%f (ang).\n", alphaVal, rcut_);
1267 chrisfen 2408 painCave.isFatal = 0;
1268     simError();
1269 chrisfen 3072 } else {
1270     alphaVal = simParams_->getDampingAlpha();
1271 chrisfen 2408 }
1272 chrisfen 3072
1273 chrisfen 2415 } else {
1274     // throw error
1275     sprintf( painCave.errMsg,
1276 gezelter 2463 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
1277     "\t(Input file specified %s .)\n"
1278     "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
1279     "or \"damped\".\n", myScreen.c_str() );
1280 chrisfen 2415 painCave.isFatal = 1;
1281     simError();
1282 chrisfen 2408 }
1283     }
1284     }
1285 chrisfen 2415
1286 chrisfen 2309 // let's pass some summation method variables to fortran
1287 chrisfen 2552 setElectrostaticSummationMethod( &esm );
1288 gezelter 2508 setFortranElectrostaticMethod( &esm );
1289 chrisfen 2408 setScreeningMethod( &sm );
1290     setDampingAlpha( &alphaVal );
1291 chrisfen 2309 setReactionFieldDielectric( &dielectric );
1292 gezelter 2463 initFortranFF( &errorOut );
1293 chrisfen 2297 }
1294    
1295 chrisfen 2425 void SimInfo::setupSwitchingFunction() {
1296     int ft = CUBIC;
1297    
1298     if (simParams_->haveSwitchingFunctionType()) {
1299     std::string funcType = simParams_->getSwitchingFunctionType();
1300     toUpper(funcType);
1301     if (funcType == "CUBIC") {
1302     ft = CUBIC;
1303     } else {
1304     if (funcType == "FIFTH_ORDER_POLYNOMIAL") {
1305     ft = FIFTH_ORDER_POLY;
1306     } else {
1307     // throw error
1308     sprintf( painCave.errMsg,
1309     "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() );
1310     painCave.isFatal = 1;
1311     simError();
1312     }
1313     }
1314     }
1315    
1316     // send switching function notification to switcheroo
1317     setFunctionType(&ft);
1318    
1319     }
1320    
1321 chrisfen 2917 void SimInfo::setupAccumulateBoxDipole() {
1322    
1323     // we only call setAccumulateBoxDipole if the accumulateBoxDipole parameter is true
1324     if ( simParams_->haveAccumulateBoxDipole() )
1325     if ( simParams_->getAccumulateBoxDipole() ) {
1326     setAccumulateBoxDipole();
1327     calcBoxDipole_ = true;
1328     }
1329    
1330     }
1331    
1332 gezelter 2204 void SimInfo::addProperty(GenericData* genData) {
1333 gezelter 1930 properties_.addProperty(genData);
1334 gezelter 2204 }
1335 gezelter 1490
1336 gezelter 2204 void SimInfo::removeProperty(const std::string& propName) {
1337 gezelter 1930 properties_.removeProperty(propName);
1338 gezelter 2204 }
1339 gezelter 1490
1340 gezelter 2204 void SimInfo::clearProperties() {
1341 gezelter 1930 properties_.clearProperties();
1342 gezelter 2204 }
1343 gezelter 1490
1344 gezelter 2204 std::vector<std::string> SimInfo::getPropertyNames() {
1345 gezelter 1930 return properties_.getPropertyNames();
1346 gezelter 2204 }
1347 gezelter 1930
1348 gezelter 2204 std::vector<GenericData*> SimInfo::getProperties() {
1349 gezelter 1930 return properties_.getProperties();
1350 gezelter 2204 }
1351 gezelter 1490
1352 gezelter 2204 GenericData* SimInfo::getPropertyByName(const std::string& propName) {
1353 gezelter 1930 return properties_.getPropertyByName(propName);
1354 gezelter 2204 }
1355 gezelter 1490
1356 gezelter 2204 void SimInfo::setSnapshotManager(SnapshotManager* sman) {
1357 tim 2116 if (sman_ == sman) {
1358 gezelter 2204 return;
1359 tim 2116 }
1360     delete sman_;
1361 gezelter 1930 sman_ = sman;
1362 gezelter 1490
1363 gezelter 1930 Molecule* mol;
1364     RigidBody* rb;
1365     Atom* atom;
1366     SimInfo::MoleculeIterator mi;
1367     Molecule::RigidBodyIterator rbIter;
1368     Molecule::AtomIterator atomIter;;
1369    
1370     for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) {
1371    
1372 gezelter 2204 for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) {
1373     atom->setSnapshotManager(sman_);
1374     }
1375 gezelter 1930
1376 gezelter 2204 for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) {
1377     rb->setSnapshotManager(sman_);
1378     }
1379 gezelter 1930 }
1380 gezelter 1490
1381 gezelter 2204 }
1382 gezelter 1490
1383 gezelter 2204 Vector3d SimInfo::getComVel(){
1384 gezelter 1930 SimInfo::MoleculeIterator i;
1385     Molecule* mol;
1386 gezelter 1490
1387 gezelter 1930 Vector3d comVel(0.0);
1388 tim 2759 RealType totalMass = 0.0;
1389 gezelter 1490
1390 gezelter 1930
1391     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1392 tim 2759 RealType mass = mol->getMass();
1393 gezelter 2204 totalMass += mass;
1394     comVel += mass * mol->getComVel();
1395 gezelter 1930 }
1396 gezelter 1490
1397 gezelter 1930 #ifdef IS_MPI
1398 tim 2759 RealType tmpMass = totalMass;
1399 gezelter 1930 Vector3d tmpComVel(comVel);
1400 tim 2759 MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1401     MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1402 gezelter 1930 #endif
1403    
1404     comVel /= totalMass;
1405    
1406     return comVel;
1407 gezelter 2204 }
1408 gezelter 1490
1409 gezelter 2204 Vector3d SimInfo::getCom(){
1410 gezelter 1930 SimInfo::MoleculeIterator i;
1411     Molecule* mol;
1412 gezelter 1490
1413 gezelter 1930 Vector3d com(0.0);
1414 tim 2759 RealType totalMass = 0.0;
1415 gezelter 1930
1416     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1417 tim 2759 RealType mass = mol->getMass();
1418 gezelter 2204 totalMass += mass;
1419     com += mass * mol->getCom();
1420 gezelter 1930 }
1421 gezelter 1490
1422     #ifdef IS_MPI
1423 tim 2759 RealType tmpMass = totalMass;
1424 gezelter 1930 Vector3d tmpCom(com);
1425 tim 2759 MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1426     MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1427 gezelter 1490 #endif
1428    
1429 gezelter 1930 com /= totalMass;
1430 gezelter 1490
1431 gezelter 1930 return com;
1432 gezelter 1490
1433 gezelter 2204 }
1434 gezelter 1930
1435 gezelter 2204 std::ostream& operator <<(std::ostream& o, SimInfo& info) {
1436 gezelter 1930
1437     return o;
1438 gezelter 2204 }
1439 chuckv 2252
1440    
1441     /*
1442     Returns center of mass and center of mass velocity in one function call.
1443     */
1444    
1445     void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){
1446     SimInfo::MoleculeIterator i;
1447     Molecule* mol;
1448    
1449    
1450 tim 2759 RealType totalMass = 0.0;
1451 chuckv 2252
1452 gezelter 1930
1453 chuckv 2252 for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1454 tim 2759 RealType mass = mol->getMass();
1455 chuckv 2252 totalMass += mass;
1456     com += mass * mol->getCom();
1457     comVel += mass * mol->getComVel();
1458     }
1459    
1460     #ifdef IS_MPI
1461 tim 2759 RealType tmpMass = totalMass;
1462 chuckv 2252 Vector3d tmpCom(com);
1463     Vector3d tmpComVel(comVel);
1464 tim 2759 MPI_Allreduce(&tmpMass,&totalMass,1,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1465     MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1466     MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1467 chuckv 2252 #endif
1468    
1469     com /= totalMass;
1470     comVel /= totalMass;
1471     }
1472    
1473     /*
1474     Return intertia tensor for entire system and angular momentum Vector.
1475 chuckv 2256
1476    
1477     [ Ixx -Ixy -Ixz ]
1478     J =| -Iyx Iyy -Iyz |
1479     [ -Izx -Iyz Izz ]
1480 chuckv 2252 */
1481    
1482     void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){
1483    
1484    
1485 tim 2759 RealType xx = 0.0;
1486     RealType yy = 0.0;
1487     RealType zz = 0.0;
1488     RealType xy = 0.0;
1489     RealType xz = 0.0;
1490     RealType yz = 0.0;
1491 chuckv 2252 Vector3d com(0.0);
1492     Vector3d comVel(0.0);
1493    
1494     getComAll(com, comVel);
1495    
1496     SimInfo::MoleculeIterator i;
1497     Molecule* mol;
1498    
1499     Vector3d thisq(0.0);
1500     Vector3d thisv(0.0);
1501    
1502 tim 2759 RealType thisMass = 0.0;
1503 chuckv 2252
1504    
1505    
1506    
1507     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1508    
1509     thisq = mol->getCom()-com;
1510     thisv = mol->getComVel()-comVel;
1511     thisMass = mol->getMass();
1512     // Compute moment of intertia coefficients.
1513     xx += thisq[0]*thisq[0]*thisMass;
1514     yy += thisq[1]*thisq[1]*thisMass;
1515     zz += thisq[2]*thisq[2]*thisMass;
1516    
1517     // compute products of intertia
1518     xy += thisq[0]*thisq[1]*thisMass;
1519     xz += thisq[0]*thisq[2]*thisMass;
1520     yz += thisq[1]*thisq[2]*thisMass;
1521    
1522     angularMomentum += cross( thisq, thisv ) * thisMass;
1523    
1524     }
1525    
1526    
1527     inertiaTensor(0,0) = yy + zz;
1528     inertiaTensor(0,1) = -xy;
1529     inertiaTensor(0,2) = -xz;
1530     inertiaTensor(1,0) = -xy;
1531 chuckv 2256 inertiaTensor(1,1) = xx + zz;
1532 chuckv 2252 inertiaTensor(1,2) = -yz;
1533     inertiaTensor(2,0) = -xz;
1534     inertiaTensor(2,1) = -yz;
1535     inertiaTensor(2,2) = xx + yy;
1536    
1537     #ifdef IS_MPI
1538     Mat3x3d tmpI(inertiaTensor);
1539     Vector3d tmpAngMom;
1540 tim 2759 MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1541     MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1542 chuckv 2252 #endif
1543    
1544     return;
1545     }
1546    
1547     //Returns the angular momentum of the system
1548     Vector3d SimInfo::getAngularMomentum(){
1549    
1550     Vector3d com(0.0);
1551     Vector3d comVel(0.0);
1552     Vector3d angularMomentum(0.0);
1553    
1554     getComAll(com,comVel);
1555    
1556     SimInfo::MoleculeIterator i;
1557     Molecule* mol;
1558    
1559 chuckv 2256 Vector3d thisr(0.0);
1560     Vector3d thisp(0.0);
1561 chuckv 2252
1562 tim 2759 RealType thisMass;
1563 chuckv 2252
1564     for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) {
1565 chuckv 2256 thisMass = mol->getMass();
1566     thisr = mol->getCom()-com;
1567     thisp = (mol->getComVel()-comVel)*thisMass;
1568 chuckv 2252
1569 chuckv 2256 angularMomentum += cross( thisr, thisp );
1570    
1571 chuckv 2252 }
1572    
1573     #ifdef IS_MPI
1574     Vector3d tmpAngMom;
1575 tim 2759 MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_REALTYPE,MPI_SUM, MPI_COMM_WORLD);
1576 chuckv 2252 #endif
1577    
1578     return angularMomentum;
1579     }
1580    
1581 tim 2982 StuntDouble* SimInfo::getIOIndexToIntegrableObject(int index) {
1582     return IOIndexToIntegrableObject.at(index);
1583     }
1584    
1585     void SimInfo::setIOIndexToIntegrableObject(const std::vector<StuntDouble*>& v) {
1586     IOIndexToIntegrableObject= v;
1587     }
1588    
1589 chuckv 3100 /* Returns the Volume of the simulation based on a ellipsoid with semi-axes
1590     based on the radius of gyration V=4/3*Pi*R_1*R_2*R_3
1591     where R_i are related to the principle inertia moments R_i = sqrt(C*I_i/N), this reduces to
1592     V = 4/3*Pi*(C/N)^3/2*sqrt(det(I)). See S.E. Baltazar et. al. Comp. Mat. Sci. 37 (2006) 526-536.
1593     */
1594     void SimInfo::getGyrationalVolume(RealType &volume){
1595     Mat3x3d intTensor;
1596     RealType det;
1597     Vector3d dummyAngMom;
1598     RealType sysconstants;
1599     RealType geomCnst;
1600    
1601     geomCnst = 3.0/2.0;
1602     /* Get the inertial tensor and angular momentum for free*/
1603     getInertiaTensor(intTensor,dummyAngMom);
1604    
1605     det = intTensor.determinant();
1606     sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1607     volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(det);
1608     return;
1609     }
1610    
1611     void SimInfo::getGyrationalVolume(RealType &volume, RealType &detI){
1612     Mat3x3d intTensor;
1613     Vector3d dummyAngMom;
1614     RealType sysconstants;
1615     RealType geomCnst;
1616    
1617     geomCnst = 3.0/2.0;
1618     /* Get the inertial tensor and angular momentum for free*/
1619     getInertiaTensor(intTensor,dummyAngMom);
1620    
1621     detI = intTensor.determinant();
1622     sysconstants = geomCnst/(RealType)nGlobalIntegrableObjects_;
1623     volume = 4.0/3.0*NumericConstant::PI*pow(sysconstants,3.0/2.0)*sqrt(detI);
1624     return;
1625     }
1626 tim 2982 /*
1627     void SimInfo::setStuntDoubleFromGlobalIndex(std::vector<StuntDouble*> v) {
1628     assert( v.size() == nAtoms_ + nRigidBodies_);
1629     sdByGlobalIndex_ = v;
1630     }
1631    
1632     StuntDouble* SimInfo::getStuntDoubleFromGlobalIndex(int index) {
1633     //assert(index < nAtoms_ + nRigidBodies_);
1634     return sdByGlobalIndex_.at(index);
1635     }
1636     */
1637 gezelter 1930 }//end namespace oopse
1638