# | Line 1 | Line 1 | |
---|---|---|
1 | < | #include <stdlib.h> |
2 | < | #include <string.h> |
3 | < | #include <math.h> |
1 | > | /* |
2 | > | * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 | > | * |
4 | > | * The University of Notre Dame grants you ("Licensee") a |
5 | > | * non-exclusive, royalty free, license to use, modify and |
6 | > | * redistribute this software in source and binary code form, provided |
7 | > | * that the following conditions are met: |
8 | > | * |
9 | > | * 1. Acknowledgement of the program authors must be made in any |
10 | > | * publication of scientific results based in part on use of the |
11 | > | * program. An acceptable form of acknowledgement is citation of |
12 | > | * the article in which the program was described (Matthew |
13 | > | * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 | > | * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 | > | * Parallel Simulation Engine for Molecular Dynamics," |
16 | > | * J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 | > | * |
18 | > | * 2. Redistributions of source code must retain the above copyright |
19 | > | * notice, this list of conditions and the following disclaimer. |
20 | > | * |
21 | > | * 3. Redistributions in binary form must reproduce the above copyright |
22 | > | * notice, this list of conditions and the following disclaimer in the |
23 | > | * documentation and/or other materials provided with the |
24 | > | * distribution. |
25 | > | * |
26 | > | * This software is provided "AS IS," without a warranty of any |
27 | > | * kind. All express or implied conditions, representations and |
28 | > | * warranties, including any implied warranty of merchantability, |
29 | > | * fitness for a particular purpose or non-infringement, are hereby |
30 | > | * excluded. The University of Notre Dame and its licensors shall not |
31 | > | * be liable for any damages suffered by licensee as a result of |
32 | > | * using, modifying or distributing the software or its |
33 | > | * derivatives. In no event will the University of Notre Dame or its |
34 | > | * licensors be liable for any lost revenue, profit or data, or for |
35 | > | * direct, indirect, special, consequential, incidental or punitive |
36 | > | * damages, however caused and regardless of the theory of liability, |
37 | > | * arising out of the use of or inability to use software, even if the |
38 | > | * University of Notre Dame has been advised of the possibility of |
39 | > | * such damages. |
40 | > | */ |
41 | > | |
42 | > | /** |
43 | > | * @file SimInfo.cpp |
44 | > | * @author tlin |
45 | > | * @date 11/02/2004 |
46 | > | * @version 1.0 |
47 | > | */ |
48 | ||
49 | < | #include <iostream> |
50 | < | using namespace std; |
49 | > | #include <algorithm> |
50 | > | #include <set> |
51 | ||
52 | #include "brains/SimInfo.hpp" | |
53 | < | #define __C |
54 | < | #include "brains/fSimulation.h" |
55 | < | #include "utils/simError.h" |
56 | < | #include "UseTheForce/DarkSide/simulation_interface.h" |
53 | > | #include "math/Vector3.hpp" |
54 | > | #include "primitives/Molecule.hpp" |
55 | > | #include "UseTheForce/fCutoffPolicy.h" |
56 | > | #include "UseTheForce/DarkSide/fElectrostaticSummationMethod.h" |
57 | > | #include "UseTheForce/DarkSide/fElectrostaticScreeningMethod.h" |
58 | > | #include "UseTheForce/DarkSide/fSwitchingFunctionType.h" |
59 | > | #include "UseTheForce/doForces_interface.h" |
60 | > | #include "UseTheForce/DarkSide/electrostatic_interface.h" |
61 | #include "UseTheForce/notifyCutoffs_interface.h" | |
62 | + | #include "UseTheForce/DarkSide/switcheroo_interface.h" |
63 | + | #include "utils/MemoryUtils.hpp" |
64 | + | #include "utils/simError.h" |
65 | + | #include "selection/SelectionManager.hpp" |
66 | ||
15 | – | //#include "UseTheForce/fortranWrappers.hpp" |
16 | – | |
17 | – | #include "math/MatVec3.h" |
18 | – | |
67 | #ifdef IS_MPI | |
68 | < | #include "brains/mpiSimulation.hpp" |
69 | < | #endif |
68 | > | #include "UseTheForce/mpiComponentPlan.h" |
69 | > | #include "UseTheForce/DarkSide/simParallel_interface.h" |
70 | > | #endif |
71 | ||
72 | < | inline double roundMe( double x ){ |
24 | < | return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); |
25 | < | } |
26 | < | |
27 | < | inline double min( double a, double b ){ |
28 | < | return (a < b ) ? a : b; |
29 | < | } |
72 | > | namespace oopse { |
73 | ||
74 | < | SimInfo* currentInfo; |
74 | > | SimInfo::SimInfo(MakeStamps* stamps, std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
75 | > | ForceField* ff, Globals* simParams) : |
76 | > | stamps_(stamps), forceField_(ff), simParams_(simParams), |
77 | > | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
78 | > | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), |
79 | > | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), |
80 | > | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), |
81 | > | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), |
82 | > | sman_(NULL), fortranInitialized_(false) { |
83 | ||
84 | < | SimInfo::SimInfo(){ |
84 | > | |
85 | > | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
86 | > | MoleculeStamp* molStamp; |
87 | > | int nMolWithSameStamp; |
88 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
89 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
90 | > | CutoffGroupStamp* cgStamp; |
91 | > | RigidBodyStamp* rbStamp; |
92 | > | int nRigidAtoms = 0; |
93 | > | |
94 | > | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
95 | > | molStamp = i->first; |
96 | > | nMolWithSameStamp = i->second; |
97 | > | |
98 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); |
99 | ||
100 | < | n_constraints = 0; |
101 | < | nZconstraints = 0; |
37 | < | n_oriented = 0; |
38 | < | n_dipoles = 0; |
39 | < | ndf = 0; |
40 | < | ndfRaw = 0; |
41 | < | nZconstraints = 0; |
42 | < | the_integrator = NULL; |
43 | < | setTemp = 0; |
44 | < | thermalTime = 0.0; |
45 | < | currentTime = 0.0; |
46 | < | rCut = 0.0; |
47 | < | rSw = 0.0; |
100 | > | //calculate atoms in molecules |
101 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
102 | ||
49 | – | haveRcut = 0; |
50 | – | haveRsw = 0; |
51 | – | boxIsInit = 0; |
52 | – | |
53 | – | resetTime = 1e99; |
103 | ||
104 | < | orthoRhombic = 0; |
105 | < | orthoTolerance = 1E-6; |
106 | < | useInitXSstate = true; |
104 | > | //calculate atoms in cutoff groups |
105 | > | int nAtomsInGroups = 0; |
106 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
107 | > | |
108 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
109 | > | cgStamp = molStamp->getCutoffGroup(j); |
110 | > | nAtomsInGroups += cgStamp->getNMembers(); |
111 | > | } |
112 | ||
113 | < | usePBC = 0; |
60 | < | useDirectionalAtoms = 0; |
61 | < | useLennardJones = 0; |
62 | < | useElectrostatics = 0; |
63 | < | useCharges = 0; |
64 | < | useDipoles = 0; |
65 | < | useSticky = 0; |
66 | < | useGayBerne = 0; |
67 | < | useEAM = 0; |
68 | < | useShapes = 0; |
69 | < | useFLARB = 0; |
113 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
114 | ||
115 | < | useSolidThermInt = 0; |
72 | < | useLiquidThermInt = 0; |
115 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
116 | ||
117 | < | haveCutoffGroups = false; |
117 | > | //calculate atoms in rigid bodies |
118 | > | int nAtomsInRigidBodies = 0; |
119 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
120 | > | |
121 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { |
122 | > | rbStamp = molStamp->getRigidBody(j); |
123 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); |
124 | > | } |
125 | ||
126 | < | excludes = Exclude::Instance(); |
126 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
127 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
128 | > | |
129 | > | } |
130 | ||
131 | < | myConfiguration = new SimState(); |
131 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff |
132 | > | //group therefore the total number of cutoff groups in the system is |
133 | > | //equal to the total number of atoms minus number of atoms belong to |
134 | > | //cutoff group defined in meta-data file plus the number of cutoff |
135 | > | //groups defined in meta-data file |
136 | > | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
137 | ||
138 | < | has_minimizer = false; |
139 | < | the_minimizer =NULL; |
138 | > | //every free atom (atom does not belong to rigid bodies) is an |
139 | > | //integrable object therefore the total number of integrable objects |
140 | > | //in the system is equal to the total number of atoms minus number of |
141 | > | //atoms belong to rigid body defined in meta-data file plus the number |
142 | > | //of rigid bodies defined in meta-data file |
143 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms |
144 | > | + nGlobalRigidBodies_; |
145 | > | |
146 | > | nGlobalMols_ = molStampIds_.size(); |
147 | ||
148 | < | ngroup = 0; |
148 | > | #ifdef IS_MPI |
149 | > | molToProcMap_.resize(nGlobalMols_); |
150 | > | #endif |
151 | ||
152 | < | } |
152 | > | } |
153 | ||
154 | + | SimInfo::~SimInfo() { |
155 | + | std::map<int, Molecule*>::iterator i; |
156 | + | for (i = molecules_.begin(); i != molecules_.end(); ++i) { |
157 | + | delete i->second; |
158 | + | } |
159 | + | molecules_.clear(); |
160 | + | |
161 | + | delete stamps_; |
162 | + | delete sman_; |
163 | + | delete simParams_; |
164 | + | delete forceField_; |
165 | + | } |
166 | ||
167 | < | SimInfo::~SimInfo(){ |
167 | > | int SimInfo::getNGlobalConstraints() { |
168 | > | int nGlobalConstraints; |
169 | > | #ifdef IS_MPI |
170 | > | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
171 | > | MPI_COMM_WORLD); |
172 | > | #else |
173 | > | nGlobalConstraints = nConstraints_; |
174 | > | #endif |
175 | > | return nGlobalConstraints; |
176 | > | } |
177 | ||
178 | < | delete myConfiguration; |
178 | > | bool SimInfo::addMolecule(Molecule* mol) { |
179 | > | MoleculeIterator i; |
180 | ||
181 | < | map<string, GenericData*>::iterator i; |
182 | < | |
94 | < | for(i = properties.begin(); i != properties.end(); i++) |
95 | < | delete (*i).second; |
181 | > | i = molecules_.find(mol->getGlobalIndex()); |
182 | > | if (i == molecules_.end() ) { |
183 | ||
184 | < | } |
184 | > | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); |
185 | > | |
186 | > | nAtoms_ += mol->getNAtoms(); |
187 | > | nBonds_ += mol->getNBonds(); |
188 | > | nBends_ += mol->getNBends(); |
189 | > | nTorsions_ += mol->getNTorsions(); |
190 | > | nRigidBodies_ += mol->getNRigidBodies(); |
191 | > | nIntegrableObjects_ += mol->getNIntegrableObjects(); |
192 | > | nCutoffGroups_ += mol->getNCutoffGroups(); |
193 | > | nConstraints_ += mol->getNConstraintPairs(); |
194 | ||
195 | < | void SimInfo::setBox(double newBox[3]) { |
196 | < | |
197 | < | int i, j; |
198 | < | double tempMat[3][3]; |
199 | < | |
104 | < | for(i=0; i<3; i++) |
105 | < | for (j=0; j<3; j++) tempMat[i][j] = 0.0;; |
106 | < | |
107 | < | tempMat[0][0] = newBox[0]; |
108 | < | tempMat[1][1] = newBox[1]; |
109 | < | tempMat[2][2] = newBox[2]; |
110 | < | |
111 | < | setBoxM( tempMat ); |
112 | < | |
113 | < | } |
114 | < | |
115 | < | void SimInfo::setBoxM( double theBox[3][3] ){ |
116 | < | |
117 | < | int i, j; |
118 | < | double FortranHmat[9]; // to preserve compatibility with Fortran the |
119 | < | // ordering in the array is as follows: |
120 | < | // [ 0 3 6 ] |
121 | < | // [ 1 4 7 ] |
122 | < | // [ 2 5 8 ] |
123 | < | double FortranHmatInv[9]; // the inverted Hmat (for Fortran); |
124 | < | |
125 | < | if( !boxIsInit ) boxIsInit = 1; |
126 | < | |
127 | < | for(i=0; i < 3; i++) |
128 | < | for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; |
129 | < | |
130 | < | calcBoxL(); |
131 | < | calcHmatInv(); |
132 | < | |
133 | < | for(i=0; i < 3; i++) { |
134 | < | for (j=0; j < 3; j++) { |
135 | < | FortranHmat[3*j + i] = Hmat[i][j]; |
136 | < | FortranHmatInv[3*j + i] = HmatInv[i][j]; |
195 | > | addExcludePairs(mol); |
196 | > | |
197 | > | return true; |
198 | > | } else { |
199 | > | return false; |
200 | } | |
201 | } | |
202 | ||
203 | < | setFortranBox(FortranHmat, FortranHmatInv, &orthoRhombic); |
204 | < | |
205 | < | } |
143 | < | |
203 | > | bool SimInfo::removeMolecule(Molecule* mol) { |
204 | > | MoleculeIterator i; |
205 | > | i = molecules_.find(mol->getGlobalIndex()); |
206 | ||
207 | < | void SimInfo::getBoxM (double theBox[3][3]) { |
207 | > | if (i != molecules_.end() ) { |
208 | ||
209 | < | int i, j; |
210 | < | for(i=0; i<3; i++) |
211 | < | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; |
212 | < | } |
209 | > | assert(mol == i->second); |
210 | > | |
211 | > | nAtoms_ -= mol->getNAtoms(); |
212 | > | nBonds_ -= mol->getNBonds(); |
213 | > | nBends_ -= mol->getNBends(); |
214 | > | nTorsions_ -= mol->getNTorsions(); |
215 | > | nRigidBodies_ -= mol->getNRigidBodies(); |
216 | > | nIntegrableObjects_ -= mol->getNIntegrableObjects(); |
217 | > | nCutoffGroups_ -= mol->getNCutoffGroups(); |
218 | > | nConstraints_ -= mol->getNConstraintPairs(); |
219 | ||
220 | + | removeExcludePairs(mol); |
221 | + | molecules_.erase(mol->getGlobalIndex()); |
222 | ||
223 | < | void SimInfo::scaleBox(double scale) { |
224 | < | double theBox[3][3]; |
225 | < | int i, j; |
223 | > | delete mol; |
224 | > | |
225 | > | return true; |
226 | > | } else { |
227 | > | return false; |
228 | > | } |
229 | ||
157 | – | // cerr << "Scaling box by " << scale << "\n"; |
230 | ||
231 | < | for(i=0; i<3; i++) |
160 | < | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; |
231 | > | } |
232 | ||
233 | < | setBoxM(theBox); |
233 | > | |
234 | > | Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { |
235 | > | i = molecules_.begin(); |
236 | > | return i == molecules_.end() ? NULL : i->second; |
237 | > | } |
238 | ||
239 | < | } |
239 | > | Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { |
240 | > | ++i; |
241 | > | return i == molecules_.end() ? NULL : i->second; |
242 | > | } |
243 | ||
166 | – | void SimInfo::calcHmatInv( void ) { |
167 | – | |
168 | – | int oldOrtho; |
169 | – | int i,j; |
170 | – | double smallDiag; |
171 | – | double tol; |
172 | – | double sanity[3][3]; |
244 | ||
245 | < | invertMat3( Hmat, HmatInv ); |
245 | > | void SimInfo::calcNdf() { |
246 | > | int ndf_local; |
247 | > | MoleculeIterator i; |
248 | > | std::vector<StuntDouble*>::iterator j; |
249 | > | Molecule* mol; |
250 | > | StuntDouble* integrableObject; |
251 | ||
252 | < | // check to see if Hmat is orthorhombic |
253 | < | |
254 | < | oldOrtho = orthoRhombic; |
252 | > | ndf_local = 0; |
253 | > | |
254 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
255 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
256 | > | integrableObject = mol->nextIntegrableObject(j)) { |
257 | ||
258 | < | smallDiag = fabs(Hmat[0][0]); |
181 | < | if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); |
182 | < | if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); |
183 | < | tol = smallDiag * orthoTolerance; |
258 | > | ndf_local += 3; |
259 | ||
260 | < | orthoRhombic = 1; |
261 | < | |
262 | < | for (i = 0; i < 3; i++ ) { |
263 | < | for (j = 0 ; j < 3; j++) { |
264 | < | if (i != j) { |
265 | < | if (orthoRhombic) { |
266 | < | if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; |
267 | < | } |
268 | < | } |
269 | < | } |
195 | < | } |
196 | < | |
197 | < | if( oldOrtho != orthoRhombic ){ |
260 | > | if (integrableObject->isDirectional()) { |
261 | > | if (integrableObject->isLinear()) { |
262 | > | ndf_local += 2; |
263 | > | } else { |
264 | > | ndf_local += 3; |
265 | > | } |
266 | > | } |
267 | > | |
268 | > | }//end for (integrableObject) |
269 | > | }// end for (mol) |
270 | ||
271 | < | if( orthoRhombic ) { |
272 | < | sprintf( painCave.errMsg, |
201 | < | "OOPSE is switching from the default Non-Orthorhombic\n" |
202 | < | "\tto the faster Orthorhombic periodic boundary computations.\n" |
203 | < | "\tThis is usually a good thing, but if you wan't the\n" |
204 | < | "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" |
205 | < | "\tvariable ( currently set to %G ) smaller.\n", |
206 | < | orthoTolerance); |
207 | < | painCave.severity = OOPSE_INFO; |
208 | < | simError(); |
209 | < | } |
210 | < | else { |
211 | < | sprintf( painCave.errMsg, |
212 | < | "OOPSE is switching from the faster Orthorhombic to the more\n" |
213 | < | "\tflexible Non-Orthorhombic periodic boundary computations.\n" |
214 | < | "\tThis is usually because the box has deformed under\n" |
215 | < | "\tNPTf integration. If you wan't to live on the edge with\n" |
216 | < | "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" |
217 | < | "\tvariable ( currently set to %G ) larger.\n", |
218 | < | orthoTolerance); |
219 | < | painCave.severity = OOPSE_WARNING; |
220 | < | simError(); |
221 | < | } |
222 | < | } |
223 | < | } |
271 | > | // n_constraints is local, so subtract them on each processor |
272 | > | ndf_local -= nConstraints_; |
273 | ||
274 | < | void SimInfo::calcBoxL( void ){ |
274 | > | #ifdef IS_MPI |
275 | > | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
276 | > | #else |
277 | > | ndf_ = ndf_local; |
278 | > | #endif |
279 | ||
280 | < | double dx, dy, dz, dsq; |
280 | > | // nZconstraints_ is global, as are the 3 COM translations for the |
281 | > | // entire system: |
282 | > | ndf_ = ndf_ - 3 - nZconstraint_; |
283 | ||
284 | < | // boxVol = Determinant of Hmat |
284 | > | } |
285 | ||
286 | < | boxVol = matDet3( Hmat ); |
286 | > | void SimInfo::calcNdfRaw() { |
287 | > | int ndfRaw_local; |
288 | ||
289 | < | // boxLx |
290 | < | |
291 | < | dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; |
292 | < | dsq = dx*dx + dy*dy + dz*dz; |
237 | < | boxL[0] = sqrt( dsq ); |
238 | < | //maxCutoff = 0.5 * boxL[0]; |
289 | > | MoleculeIterator i; |
290 | > | std::vector<StuntDouble*>::iterator j; |
291 | > | Molecule* mol; |
292 | > | StuntDouble* integrableObject; |
293 | ||
294 | < | // boxLy |
295 | < | |
296 | < | dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; |
297 | < | dsq = dx*dx + dy*dy + dz*dz; |
298 | < | boxL[1] = sqrt( dsq ); |
299 | < | //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1]; |
294 | > | // Raw degrees of freedom that we have to set |
295 | > | ndfRaw_local = 0; |
296 | > | |
297 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
298 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
299 | > | integrableObject = mol->nextIntegrableObject(j)) { |
300 | ||
301 | < | |
248 | < | // boxLz |
249 | < | |
250 | < | dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; |
251 | < | dsq = dx*dx + dy*dy + dz*dz; |
252 | < | boxL[2] = sqrt( dsq ); |
253 | < | //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2]; |
254 | < | |
255 | < | //calculate the max cutoff |
256 | < | maxCutoff = calcMaxCutOff(); |
257 | < | |
258 | < | checkCutOffs(); |
259 | < | |
260 | < | } |
261 | < | |
262 | < | |
263 | < | double SimInfo::calcMaxCutOff(){ |
301 | > | ndfRaw_local += 3; |
302 | ||
303 | < | double ri[3], rj[3], rk[3]; |
304 | < | double rij[3], rjk[3], rki[3]; |
305 | < | double minDist; |
303 | > | if (integrableObject->isDirectional()) { |
304 | > | if (integrableObject->isLinear()) { |
305 | > | ndfRaw_local += 2; |
306 | > | } else { |
307 | > | ndfRaw_local += 3; |
308 | > | } |
309 | > | } |
310 | > | |
311 | > | } |
312 | > | } |
313 | > | |
314 | > | #ifdef IS_MPI |
315 | > | MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
316 | > | #else |
317 | > | ndfRaw_ = ndfRaw_local; |
318 | > | #endif |
319 | > | } |
320 | ||
321 | < | ri[0] = Hmat[0][0]; |
322 | < | ri[1] = Hmat[1][0]; |
271 | < | ri[2] = Hmat[2][0]; |
321 | > | void SimInfo::calcNdfTrans() { |
322 | > | int ndfTrans_local; |
323 | ||
324 | < | rj[0] = Hmat[0][1]; |
274 | < | rj[1] = Hmat[1][1]; |
275 | < | rj[2] = Hmat[2][1]; |
324 | > | ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; |
325 | ||
277 | – | rk[0] = Hmat[0][2]; |
278 | – | rk[1] = Hmat[1][2]; |
279 | – | rk[2] = Hmat[2][2]; |
280 | – | |
281 | – | crossProduct3(ri, rj, rij); |
282 | – | distXY = dotProduct3(rk,rij) / norm3(rij); |
326 | ||
327 | < | crossProduct3(rj,rk, rjk); |
328 | < | distYZ = dotProduct3(ri,rjk) / norm3(rjk); |
327 | > | #ifdef IS_MPI |
328 | > | MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
329 | > | #else |
330 | > | ndfTrans_ = ndfTrans_local; |
331 | > | #endif |
332 | ||
333 | < | crossProduct3(rk,ri, rki); |
334 | < | distZX = dotProduct3(rj,rki) / norm3(rki); |
333 | > | ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; |
334 | > | |
335 | > | } |
336 | ||
337 | < | minDist = min(min(distXY, distYZ), distZX); |
338 | < | return minDist/2; |
339 | < | |
340 | < | } |
337 | > | void SimInfo::addExcludePairs(Molecule* mol) { |
338 | > | std::vector<Bond*>::iterator bondIter; |
339 | > | std::vector<Bend*>::iterator bendIter; |
340 | > | std::vector<Torsion*>::iterator torsionIter; |
341 | > | Bond* bond; |
342 | > | Bend* bend; |
343 | > | Torsion* torsion; |
344 | > | int a; |
345 | > | int b; |
346 | > | int c; |
347 | > | int d; |
348 | > | |
349 | > | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
350 | > | a = bond->getAtomA()->getGlobalIndex(); |
351 | > | b = bond->getAtomB()->getGlobalIndex(); |
352 | > | exclude_.addPair(a, b); |
353 | > | } |
354 | ||
355 | < | void SimInfo::wrapVector( double thePos[3] ){ |
355 | > | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
356 | > | a = bend->getAtomA()->getGlobalIndex(); |
357 | > | b = bend->getAtomB()->getGlobalIndex(); |
358 | > | c = bend->getAtomC()->getGlobalIndex(); |
359 | ||
360 | < | int i; |
361 | < | double scaled[3]; |
360 | > | exclude_.addPair(a, b); |
361 | > | exclude_.addPair(a, c); |
362 | > | exclude_.addPair(b, c); |
363 | > | } |
364 | ||
365 | < | if( !orthoRhombic ){ |
366 | < | // calc the scaled coordinates. |
367 | < | |
365 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
366 | > | a = torsion->getAtomA()->getGlobalIndex(); |
367 | > | b = torsion->getAtomB()->getGlobalIndex(); |
368 | > | c = torsion->getAtomC()->getGlobalIndex(); |
369 | > | d = torsion->getAtomD()->getGlobalIndex(); |
370 | ||
371 | < | matVecMul3(HmatInv, thePos, scaled); |
372 | < | |
373 | < | for(i=0; i<3; i++) |
374 | < | scaled[i] -= roundMe(scaled[i]); |
375 | < | |
376 | < | // calc the wrapped real coordinates from the wrapped scaled coordinates |
377 | < | |
311 | < | matVecMul3(Hmat, scaled, thePos); |
371 | > | exclude_.addPair(a, b); |
372 | > | exclude_.addPair(a, c); |
373 | > | exclude_.addPair(a, d); |
374 | > | exclude_.addPair(b, c); |
375 | > | exclude_.addPair(b, d); |
376 | > | exclude_.addPair(c, d); |
377 | > | } |
378 | ||
379 | + | Molecule::RigidBodyIterator rbIter; |
380 | + | RigidBody* rb; |
381 | + | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
382 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
383 | + | for (int i = 0; i < atoms.size() -1 ; ++i) { |
384 | + | for (int j = i + 1; j < atoms.size(); ++j) { |
385 | + | a = atoms[i]->getGlobalIndex(); |
386 | + | b = atoms[j]->getGlobalIndex(); |
387 | + | exclude_.addPair(a, b); |
388 | + | } |
389 | + | } |
390 | + | } |
391 | + | |
392 | } | |
393 | < | else{ |
394 | < | // calc the scaled coordinates. |
393 | > | |
394 | > | void SimInfo::removeExcludePairs(Molecule* mol) { |
395 | > | std::vector<Bond*>::iterator bondIter; |
396 | > | std::vector<Bend*>::iterator bendIter; |
397 | > | std::vector<Torsion*>::iterator torsionIter; |
398 | > | Bond* bond; |
399 | > | Bend* bend; |
400 | > | Torsion* torsion; |
401 | > | int a; |
402 | > | int b; |
403 | > | int c; |
404 | > | int d; |
405 | ||
406 | < | for(i=0; i<3; i++) |
407 | < | scaled[i] = thePos[i]*HmatInv[i][i]; |
408 | < | |
409 | < | // wrap the scaled coordinates |
410 | < | |
322 | < | for(i=0; i<3; i++) |
323 | < | scaled[i] -= roundMe(scaled[i]); |
324 | < | |
325 | < | // calc the wrapped real coordinates from the wrapped scaled coordinates |
326 | < | |
327 | < | for(i=0; i<3; i++) |
328 | < | thePos[i] = scaled[i]*Hmat[i][i]; |
329 | < | } |
330 | < | |
331 | < | } |
406 | > | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
407 | > | a = bond->getAtomA()->getGlobalIndex(); |
408 | > | b = bond->getAtomB()->getGlobalIndex(); |
409 | > | exclude_.removePair(a, b); |
410 | > | } |
411 | ||
412 | + | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
413 | + | a = bend->getAtomA()->getGlobalIndex(); |
414 | + | b = bend->getAtomB()->getGlobalIndex(); |
415 | + | c = bend->getAtomC()->getGlobalIndex(); |
416 | ||
417 | < | int SimInfo::getNDF(){ |
418 | < | int ndf_local; |
417 | > | exclude_.removePair(a, b); |
418 | > | exclude_.removePair(a, c); |
419 | > | exclude_.removePair(b, c); |
420 | > | } |
421 | ||
422 | < | ndf_local = 0; |
423 | < | |
424 | < | for(int i = 0; i < integrableObjects.size(); i++){ |
425 | < | ndf_local += 3; |
426 | < | if (integrableObjects[i]->isDirectional()) { |
427 | < | if (integrableObjects[i]->isLinear()) |
428 | < | ndf_local += 2; |
429 | < | else |
430 | < | ndf_local += 3; |
422 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
423 | > | a = torsion->getAtomA()->getGlobalIndex(); |
424 | > | b = torsion->getAtomB()->getGlobalIndex(); |
425 | > | c = torsion->getAtomC()->getGlobalIndex(); |
426 | > | d = torsion->getAtomD()->getGlobalIndex(); |
427 | > | |
428 | > | exclude_.removePair(a, b); |
429 | > | exclude_.removePair(a, c); |
430 | > | exclude_.removePair(a, d); |
431 | > | exclude_.removePair(b, c); |
432 | > | exclude_.removePair(b, d); |
433 | > | exclude_.removePair(c, d); |
434 | } | |
435 | + | |
436 | + | Molecule::RigidBodyIterator rbIter; |
437 | + | RigidBody* rb; |
438 | + | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
439 | + | std::vector<Atom*> atoms = rb->getAtoms(); |
440 | + | for (int i = 0; i < atoms.size() -1 ; ++i) { |
441 | + | for (int j = i + 1; j < atoms.size(); ++j) { |
442 | + | a = atoms[i]->getGlobalIndex(); |
443 | + | b = atoms[j]->getGlobalIndex(); |
444 | + | exclude_.removePair(a, b); |
445 | + | } |
446 | + | } |
447 | + | } |
448 | + | |
449 | } | |
450 | ||
349 | – | // n_constraints is local, so subtract them on each processor: |
451 | ||
452 | < | ndf_local -= n_constraints; |
452 | > | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { |
453 | > | int curStampId; |
454 | ||
455 | + | //index from 0 |
456 | + | curStampId = moleculeStamps_.size(); |
457 | + | |
458 | + | moleculeStamps_.push_back(molStamp); |
459 | + | molStampIds_.insert(molStampIds_.end(), nmol, curStampId); |
460 | + | } |
461 | + | |
462 | + | void SimInfo::update() { |
463 | + | |
464 | + | setupSimType(); |
465 | + | |
466 | #ifdef IS_MPI | |
467 | < | MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
355 | < | #else |
356 | < | ndf = ndf_local; |
467 | > | setupFortranParallel(); |
468 | #endif | |
469 | ||
470 | < | // nZconstraints is global, as are the 3 COM translations for the |
360 | < | // entire system: |
470 | > | setupFortranSim(); |
471 | ||
472 | < | ndf = ndf - 3 - nZconstraints; |
472 | > | //setup fortran force field |
473 | > | /** @deprecate */ |
474 | > | int isError = 0; |
475 | > | |
476 | > | setupElectrostaticSummationMethod( isError ); |
477 | > | setupSwitchingFunction(); |
478 | ||
479 | < | return ndf; |
480 | < | } |
479 | > | if(isError){ |
480 | > | sprintf( painCave.errMsg, |
481 | > | "ForceField error: There was an error initializing the forceField in fortran.\n" ); |
482 | > | painCave.isFatal = 1; |
483 | > | simError(); |
484 | > | } |
485 | > | |
486 | > | |
487 | > | setupCutoff(); |
488 | ||
489 | < | int SimInfo::getNDFraw() { |
490 | < | int ndfRaw_local; |
489 | > | calcNdf(); |
490 | > | calcNdfRaw(); |
491 | > | calcNdfTrans(); |
492 | ||
493 | < | // Raw degrees of freedom that we have to set |
494 | < | ndfRaw_local = 0; |
493 | > | fortranInitialized_ = true; |
494 | > | } |
495 | ||
496 | < | for(int i = 0; i < integrableObjects.size(); i++){ |
497 | < | ndfRaw_local += 3; |
498 | < | if (integrableObjects[i]->isDirectional()) { |
499 | < | if (integrableObjects[i]->isLinear()) |
500 | < | ndfRaw_local += 2; |
501 | < | else |
502 | < | ndfRaw_local += 3; |
496 | > | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
497 | > | SimInfo::MoleculeIterator mi; |
498 | > | Molecule* mol; |
499 | > | Molecule::AtomIterator ai; |
500 | > | Atom* atom; |
501 | > | std::set<AtomType*> atomTypes; |
502 | > | |
503 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
504 | > | |
505 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
506 | > | atomTypes.insert(atom->getAtomType()); |
507 | > | } |
508 | > | |
509 | } | |
510 | + | |
511 | + | return atomTypes; |
512 | } | |
513 | + | |
514 | + | void SimInfo::setupSimType() { |
515 | + | std::set<AtomType*>::iterator i; |
516 | + | std::set<AtomType*> atomTypes; |
517 | + | atomTypes = getUniqueAtomTypes(); |
518 | ||
519 | < | #ifdef IS_MPI |
520 | < | MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
521 | < | #else |
522 | < | ndfRaw = ndfRaw_local; |
523 | < | #endif |
519 | > | int useLennardJones = 0; |
520 | > | int useElectrostatic = 0; |
521 | > | int useEAM = 0; |
522 | > | int useCharge = 0; |
523 | > | int useDirectional = 0; |
524 | > | int useDipole = 0; |
525 | > | int useGayBerne = 0; |
526 | > | int useSticky = 0; |
527 | > | int useStickyPower = 0; |
528 | > | int useShape = 0; |
529 | > | int useFLARB = 0; //it is not in AtomType yet |
530 | > | int useDirectionalAtom = 0; |
531 | > | int useElectrostatics = 0; |
532 | > | //usePBC and useRF are from simParams |
533 | > | int usePBC = simParams_->getUsePeriodicBoundaryConditions(); |
534 | > | int useRF; |
535 | > | int useSF; |
536 | > | std::string myMethod; |
537 | ||
538 | < | return ndfRaw; |
539 | < | } |
538 | > | // set the useRF logical |
539 | > | useRF = 0; |
540 | > | useSF = 0; |
541 | ||
392 | – | int SimInfo::getNDFtranslational() { |
393 | – | int ndfTrans_local; |
542 | ||
543 | < | ndfTrans_local = 3 * integrableObjects.size() - n_constraints; |
543 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
544 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
545 | > | toUpper(myMethod); |
546 | > | if (myMethod == "REACTION_FIELD") { |
547 | > | useRF=1; |
548 | > | } else { |
549 | > | if (myMethod == "SHIFTED_FORCE") { |
550 | > | useSF = 1; |
551 | > | } |
552 | > | } |
553 | > | } |
554 | ||
555 | + | //loop over all of the atom types |
556 | + | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
557 | + | useLennardJones |= (*i)->isLennardJones(); |
558 | + | useElectrostatic |= (*i)->isElectrostatic(); |
559 | + | useEAM |= (*i)->isEAM(); |
560 | + | useCharge |= (*i)->isCharge(); |
561 | + | useDirectional |= (*i)->isDirectional(); |
562 | + | useDipole |= (*i)->isDipole(); |
563 | + | useGayBerne |= (*i)->isGayBerne(); |
564 | + | useSticky |= (*i)->isSticky(); |
565 | + | useStickyPower |= (*i)->isStickyPower(); |
566 | + | useShape |= (*i)->isShape(); |
567 | + | } |
568 | ||
569 | < | #ifdef IS_MPI |
570 | < | MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
571 | < | #else |
401 | < | ndfTrans = ndfTrans_local; |
402 | < | #endif |
569 | > | if (useSticky || useStickyPower || useDipole || useGayBerne || useShape) { |
570 | > | useDirectionalAtom = 1; |
571 | > | } |
572 | ||
573 | < | ndfTrans = ndfTrans - 3 - nZconstraints; |
573 | > | if (useCharge || useDipole) { |
574 | > | useElectrostatics = 1; |
575 | > | } |
576 | ||
577 | < | return ndfTrans; |
578 | < | } |
577 | > | #ifdef IS_MPI |
578 | > | int temp; |
579 | ||
580 | < | int SimInfo::getTotIntegrableObjects() { |
581 | < | int nObjs_local; |
411 | < | int nObjs; |
580 | > | temp = usePBC; |
581 | > | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
582 | ||
583 | < | nObjs_local = integrableObjects.size(); |
583 | > | temp = useDirectionalAtom; |
584 | > | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
585 | ||
586 | + | temp = useLennardJones; |
587 | + | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
588 | ||
589 | < | #ifdef IS_MPI |
590 | < | MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
418 | < | #else |
419 | < | nObjs = nObjs_local; |
420 | < | #endif |
589 | > | temp = useElectrostatics; |
590 | > | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
591 | ||
592 | + | temp = useCharge; |
593 | + | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
594 | ||
595 | < | return nObjs; |
596 | < | } |
595 | > | temp = useDipole; |
596 | > | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
597 | ||
598 | < | void SimInfo::refreshSim(){ |
598 | > | temp = useSticky; |
599 | > | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
600 | ||
601 | < | simtype fInfo; |
602 | < | int isError; |
603 | < | int n_global; |
604 | < | int* excl; |
601 | > | temp = useStickyPower; |
602 | > | MPI_Allreduce(&temp, &useStickyPower, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
603 | > | |
604 | > | temp = useGayBerne; |
605 | > | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
606 | ||
607 | < | fInfo.dielect = 0.0; |
607 | > | temp = useEAM; |
608 | > | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
609 | ||
610 | < | if( useDipoles ){ |
611 | < | if( useReactionField )fInfo.dielect = dielectric; |
437 | < | } |
610 | > | temp = useShape; |
611 | > | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
612 | ||
613 | < | fInfo.SIM_uses_PBC = usePBC; |
613 | > | temp = useFLARB; |
614 | > | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
615 | ||
616 | < | if (useSticky || useDipoles || useGayBerne || useShapes) { |
617 | < | useDirectionalAtoms = 1; |
443 | < | fInfo.SIM_uses_DirectionalAtoms = useDirectionalAtoms; |
444 | < | } |
616 | > | temp = useRF; |
617 | > | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
618 | ||
619 | < | fInfo.SIM_uses_LennardJones = useLennardJones; |
619 | > | temp = useSF; |
620 | > | MPI_Allreduce(&temp, &useSF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
621 | ||
622 | < | if (useCharges || useDipoles) { |
623 | < | useElectrostatics = 1; |
624 | < | fInfo.SIM_uses_Electrostatics = useElectrostatics; |
622 | > | #endif |
623 | > | |
624 | > | fInfo_.SIM_uses_PBC = usePBC; |
625 | > | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; |
626 | > | fInfo_.SIM_uses_LennardJones = useLennardJones; |
627 | > | fInfo_.SIM_uses_Electrostatics = useElectrostatics; |
628 | > | fInfo_.SIM_uses_Charges = useCharge; |
629 | > | fInfo_.SIM_uses_Dipoles = useDipole; |
630 | > | fInfo_.SIM_uses_Sticky = useSticky; |
631 | > | fInfo_.SIM_uses_StickyPower = useStickyPower; |
632 | > | fInfo_.SIM_uses_GayBerne = useGayBerne; |
633 | > | fInfo_.SIM_uses_EAM = useEAM; |
634 | > | fInfo_.SIM_uses_Shapes = useShape; |
635 | > | fInfo_.SIM_uses_FLARB = useFLARB; |
636 | > | fInfo_.SIM_uses_RF = useRF; |
637 | > | fInfo_.SIM_uses_SF = useSF; |
638 | > | |
639 | > | if( myMethod == "REACTION_FIELD") { |
640 | > | |
641 | > | if (simParams_->haveDielectric()) { |
642 | > | fInfo_.dielect = simParams_->getDielectric(); |
643 | > | } else { |
644 | > | sprintf(painCave.errMsg, |
645 | > | "SimSetup Error: No Dielectric constant was set.\n" |
646 | > | "\tYou are trying to use Reaction Field without" |
647 | > | "\tsetting a dielectric constant!\n"); |
648 | > | painCave.isFatal = 1; |
649 | > | simError(); |
650 | > | } |
651 | > | } |
652 | > | |
653 | } | |
654 | ||
655 | < | fInfo.SIM_uses_Charges = useCharges; |
656 | < | fInfo.SIM_uses_Dipoles = useDipoles; |
657 | < | fInfo.SIM_uses_Sticky = useSticky; |
658 | < | fInfo.SIM_uses_GayBerne = useGayBerne; |
659 | < | fInfo.SIM_uses_EAM = useEAM; |
660 | < | fInfo.SIM_uses_Shapes = useShapes; |
661 | < | fInfo.SIM_uses_FLARB = useFLARB; |
460 | < | fInfo.SIM_uses_RF = useReactionField; |
655 | > | void SimInfo::setupFortranSim() { |
656 | > | int isError; |
657 | > | int nExclude; |
658 | > | std::vector<int> fortranGlobalGroupMembership; |
659 | > | |
660 | > | nExclude = exclude_.getSize(); |
661 | > | isError = 0; |
662 | ||
663 | < | n_exclude = excludes->getSize(); |
664 | < | excl = excludes->getFortranArray(); |
665 | < | |
666 | < | #ifdef IS_MPI |
466 | < | n_global = mpiSim->getNAtomsGlobal(); |
467 | < | #else |
468 | < | n_global = n_atoms; |
469 | < | #endif |
470 | < | |
471 | < | isError = 0; |
472 | < | |
473 | < | getFortranGroupArrays(this, FglobalGroupMembership, mfact); |
474 | < | //it may not be a good idea to pass the address of first element in vector |
475 | < | //since c++ standard does not require vector to be stored continuously in meomory |
476 | < | //Most of the compilers will organize the memory of vector continuously |
477 | < | setFortranSim( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, |
478 | < | &nGlobalExcludes, globalExcludes, molMembershipArray, |
479 | < | &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); |
663 | > | //globalGroupMembership_ is filled by SimCreator |
664 | > | for (int i = 0; i < nGlobalAtoms_; i++) { |
665 | > | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
666 | > | } |
667 | ||
668 | < | if( isError ){ |
668 | > | //calculate mass ratio of cutoff group |
669 | > | std::vector<double> mfact; |
670 | > | SimInfo::MoleculeIterator mi; |
671 | > | Molecule* mol; |
672 | > | Molecule::CutoffGroupIterator ci; |
673 | > | CutoffGroup* cg; |
674 | > | Molecule::AtomIterator ai; |
675 | > | Atom* atom; |
676 | > | double totalMass; |
677 | > | |
678 | > | //to avoid memory reallocation, reserve enough space for mfact |
679 | > | mfact.reserve(getNCutoffGroups()); |
680 | ||
681 | < | sprintf( painCave.errMsg, |
682 | < | "There was an error setting the simulation information in fortran.\n" ); |
485 | < | painCave.isFatal = 1; |
486 | < | painCave.severity = OOPSE_ERROR; |
487 | < | simError(); |
488 | < | } |
489 | < | |
490 | < | #ifdef IS_MPI |
491 | < | sprintf( checkPointMsg, |
492 | < | "succesfully sent the simulation information to fortran.\n"); |
493 | < | MPIcheckPoint(); |
494 | < | #endif // is_mpi |
495 | < | |
496 | < | this->ndf = this->getNDF(); |
497 | < | this->ndfRaw = this->getNDFraw(); |
498 | < | this->ndfTrans = this->getNDFtranslational(); |
499 | < | } |
681 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
682 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
683 | ||
684 | < | void SimInfo::setDefaultRcut( double theRcut ){ |
685 | < | |
686 | < | haveRcut = 1; |
687 | < | rCut = theRcut; |
688 | < | rList = rCut + 1.0; |
689 | < | |
690 | < | notifyFortranCutoffs( &rCut, &rSw, &rList ); |
691 | < | } |
684 | > | totalMass = cg->getMass(); |
685 | > | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { |
686 | > | // Check for massless groups - set mfact to 1 if true |
687 | > | if (totalMass != 0) |
688 | > | mfact.push_back(atom->getMass()/totalMass); |
689 | > | else |
690 | > | mfact.push_back( 1.0 ); |
691 | > | } |
692 | ||
693 | < | void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ |
693 | > | } |
694 | > | } |
695 | ||
696 | < | rSw = theRsw; |
697 | < | setDefaultRcut( theRcut ); |
514 | < | } |
696 | > | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
697 | > | std::vector<int> identArray; |
698 | ||
699 | + | //to avoid memory reallocation, reserve enough space identArray |
700 | + | identArray.reserve(getNAtoms()); |
701 | + | |
702 | + | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
703 | + | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
704 | + | identArray.push_back(atom->getIdent()); |
705 | + | } |
706 | + | } |
707 | ||
708 | < | void SimInfo::checkCutOffs( void ){ |
709 | < | |
710 | < | if( boxIsInit ){ |
708 | > | //fill molMembershipArray |
709 | > | //molMembershipArray is filled by SimCreator |
710 | > | std::vector<int> molMembershipArray(nGlobalAtoms_); |
711 | > | for (int i = 0; i < nGlobalAtoms_; i++) { |
712 | > | molMembershipArray[i] = globalMolMembership_[i] + 1; |
713 | > | } |
714 | ||
715 | < | //we need to check cutOffs against the box |
716 | < | |
717 | < | if( rCut > maxCutoff ){ |
715 | > | //setup fortran simulation |
716 | > | int nGlobalExcludes = 0; |
717 | > | int* globalExcludes = NULL; |
718 | > | int* excludeList = exclude_.getExcludeList(); |
719 | > | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , |
720 | > | &nGlobalExcludes, globalExcludes, &molMembershipArray[0], |
721 | > | &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); |
722 | > | |
723 | > | if( isError ){ |
724 | > | |
725 | sprintf( painCave.errMsg, | |
726 | < | "cutoffRadius is too large for the current periodic box.\n" |
526 | < | "\tCurrent Value of cutoffRadius = %G at time %G\n " |
527 | < | "\tThis is larger than half of at least one of the\n" |
528 | < | "\tperiodic box vectors. Right now, the Box matrix is:\n" |
529 | < | "\n" |
530 | < | "\t[ %G %G %G ]\n" |
531 | < | "\t[ %G %G %G ]\n" |
532 | < | "\t[ %G %G %G ]\n", |
533 | < | rCut, currentTime, |
534 | < | Hmat[0][0], Hmat[0][1], Hmat[0][2], |
535 | < | Hmat[1][0], Hmat[1][1], Hmat[1][2], |
536 | < | Hmat[2][0], Hmat[2][1], Hmat[2][2]); |
537 | < | painCave.severity = OOPSE_ERROR; |
726 | > | "There was an error setting the simulation information in fortran.\n" ); |
727 | painCave.isFatal = 1; | |
728 | + | painCave.severity = OOPSE_ERROR; |
729 | simError(); | |
730 | < | } |
731 | < | } else { |
732 | < | // initialize this stuff before using it, OK? |
733 | < | sprintf( painCave.errMsg, |
734 | < | "Trying to check cutoffs without a box.\n" |
735 | < | "\tOOPSE should have better programmers than that.\n" ); |
736 | < | painCave.severity = OOPSE_ERROR; |
547 | < | painCave.isFatal = 1; |
548 | < | simError(); |
730 | > | } |
731 | > | |
732 | > | #ifdef IS_MPI |
733 | > | sprintf( checkPointMsg, |
734 | > | "succesfully sent the simulation information to fortran.\n"); |
735 | > | MPIcheckPoint(); |
736 | > | #endif // is_mpi |
737 | } | |
550 | – | |
551 | – | } |
738 | ||
553 | – | void SimInfo::addProperty(GenericData* prop){ |
739 | ||
740 | < | map<string, GenericData*>::iterator result; |
741 | < | result = properties.find(prop->getID()); |
557 | < | |
558 | < | //we can't simply use properties[prop->getID()] = prop, |
559 | < | //it will cause memory leak if we already contain a propery which has the same name of prop |
560 | < | |
561 | < | if(result != properties.end()){ |
740 | > | #ifdef IS_MPI |
741 | > | void SimInfo::setupFortranParallel() { |
742 | ||
743 | < | delete (*result).second; |
744 | < | (*result).second = prop; |
745 | < | |
746 | < | } |
747 | < | else{ |
743 | > | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex |
744 | > | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
745 | > | std::vector<int> localToGlobalCutoffGroupIndex; |
746 | > | SimInfo::MoleculeIterator mi; |
747 | > | Molecule::AtomIterator ai; |
748 | > | Molecule::CutoffGroupIterator ci; |
749 | > | Molecule* mol; |
750 | > | Atom* atom; |
751 | > | CutoffGroup* cg; |
752 | > | mpiSimData parallelData; |
753 | > | int isError; |
754 | ||
755 | < | properties[prop->getID()] = prop; |
755 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
756 | ||
757 | + | //local index(index in DataStorge) of atom is important |
758 | + | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
759 | + | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
760 | + | } |
761 | + | |
762 | + | //local index of cutoff group is trivial, it only depends on the order of travesing |
763 | + | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
764 | + | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
765 | + | } |
766 | + | |
767 | + | } |
768 | + | |
769 | + | //fill up mpiSimData struct |
770 | + | parallelData.nMolGlobal = getNGlobalMolecules(); |
771 | + | parallelData.nMolLocal = getNMolecules(); |
772 | + | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
773 | + | parallelData.nAtomsLocal = getNAtoms(); |
774 | + | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
775 | + | parallelData.nGroupsLocal = getNCutoffGroups(); |
776 | + | parallelData.myNode = worldRank; |
777 | + | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
778 | + | |
779 | + | //pass mpiSimData struct and index arrays to fortran |
780 | + | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
781 | + | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
782 | + | &localToGlobalCutoffGroupIndex[0], &isError); |
783 | + | |
784 | + | if (isError) { |
785 | + | sprintf(painCave.errMsg, |
786 | + | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
787 | + | painCave.isFatal = 1; |
788 | + | simError(); |
789 | + | } |
790 | + | |
791 | + | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
792 | + | MPIcheckPoint(); |
793 | + | |
794 | + | |
795 | } | |
572 | – | |
573 | – | } |
796 | ||
797 | < | GenericData* SimInfo::getProperty(const string& propName){ |
576 | < | |
577 | < | map<string, GenericData*>::iterator result; |
578 | < | |
579 | < | //string lowerCaseName = (); |
580 | < | |
581 | < | result = properties.find(propName); |
582 | < | |
583 | < | if(result != properties.end()) |
584 | < | return (*result).second; |
585 | < | else |
586 | < | return NULL; |
587 | < | } |
797 | > | #endif |
798 | ||
799 | + | double SimInfo::calcMaxCutoffRadius() { |
800 | ||
590 | – | void SimInfo::getFortranGroupArrays(SimInfo* info, |
591 | – | vector<int>& FglobalGroupMembership, |
592 | – | vector<double>& mfact){ |
593 | – | |
594 | – | Molecule* myMols; |
595 | – | Atom** myAtoms; |
596 | – | int numAtom; |
597 | – | double mtot; |
598 | – | int numMol; |
599 | – | int numCutoffGroups; |
600 | – | CutoffGroup* myCutoffGroup; |
601 | – | vector<CutoffGroup*>::iterator iterCutoff; |
602 | – | Atom* cutoffAtom; |
603 | – | vector<Atom*>::iterator iterAtom; |
604 | – | int atomIndex; |
605 | – | double totalMass; |
606 | – | |
607 | – | mfact.clear(); |
608 | – | FglobalGroupMembership.clear(); |
609 | – | |
801 | ||
802 | < | // Fix the silly fortran indexing problem |
802 | > | std::set<AtomType*> atomTypes; |
803 | > | std::set<AtomType*>::iterator i; |
804 | > | std::vector<double> cutoffRadius; |
805 | > | |
806 | > | //get the unique atom types |
807 | > | atomTypes = getUniqueAtomTypes(); |
808 | > | |
809 | > | //query the max cutoff radius among these atom types |
810 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
811 | > | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
812 | > | } |
813 | > | |
814 | > | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); |
815 | #ifdef IS_MPI | |
816 | < | numAtom = mpiSim->getNAtomsGlobal(); |
614 | < | #else |
615 | < | numAtom = n_atoms; |
816 | > | //pick the max cutoff radius among the processors |
817 | #endif | |
617 | – | for (int i = 0; i < numAtom; i++) |
618 | – | FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); |
619 | – | |
818 | ||
819 | < | myMols = info->molecules; |
820 | < | numMol = info->n_mol; |
623 | < | for(int i = 0; i < numMol; i++){ |
624 | < | numCutoffGroups = myMols[i].getNCutoffGroups(); |
625 | < | for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); |
626 | < | myCutoffGroup != NULL; |
627 | < | myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ |
819 | > | return maxCutoffRadius; |
820 | > | } |
821 | ||
822 | < | totalMass = myCutoffGroup->getMass(); |
823 | < | |
824 | < | for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); |
825 | < | cutoffAtom != NULL; |
826 | < | cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ |
827 | < | mfact.push_back(cutoffAtom->getMass()/totalMass); |
828 | < | } |
822 | > | void SimInfo::getCutoff(double& rcut, double& rsw) { |
823 | > | |
824 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
825 | > | |
826 | > | if (!simParams_->haveCutoffRadius()){ |
827 | > | sprintf(painCave.errMsg, |
828 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
829 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
830 | > | "\tfor the cutoffRadius.\n"); |
831 | > | painCave.isFatal = 0; |
832 | > | simError(); |
833 | > | rcut = 15.0; |
834 | > | } else{ |
835 | > | rcut = simParams_->getCutoffRadius(); |
836 | > | } |
837 | > | |
838 | > | if (!simParams_->haveSwitchingRadius()){ |
839 | > | sprintf(painCave.errMsg, |
840 | > | "SimCreator Warning: No value was set for switchingRadius.\n" |
841 | > | "\tOOPSE will use a default value of\n" |
842 | > | "\t0.85 * cutoffRadius for the switchingRadius\n"); |
843 | > | painCave.isFatal = 0; |
844 | > | simError(); |
845 | > | rsw = 0.85 * rcut; |
846 | > | } else{ |
847 | > | rsw = simParams_->getSwitchingRadius(); |
848 | > | } |
849 | > | |
850 | > | } else { |
851 | > | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
852 | > | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
853 | > | |
854 | > | if (simParams_->haveCutoffRadius()) { |
855 | > | rcut = simParams_->getCutoffRadius(); |
856 | > | } else { |
857 | > | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
858 | > | rcut = calcMaxCutoffRadius(); |
859 | > | } |
860 | > | |
861 | > | if (simParams_->haveSwitchingRadius()) { |
862 | > | rsw = simParams_->getSwitchingRadius(); |
863 | > | } else { |
864 | > | rsw = rcut; |
865 | > | } |
866 | > | |
867 | } | |
868 | } | |
869 | ||
870 | < | } |
870 | > | void SimInfo::setupCutoff() { |
871 | > | getCutoff(rcut_, rsw_); |
872 | > | double rnblist = rcut_ + 1; // skin of neighbor list |
873 | > | |
874 | > | //Pass these cutoff radius etc. to fortran. This function should be called once and only once |
875 | > | |
876 | > | int cp = TRADITIONAL_CUTOFF_POLICY; |
877 | > | if (simParams_->haveCutoffPolicy()) { |
878 | > | std::string myPolicy = simParams_->getCutoffPolicy(); |
879 | > | toUpper(myPolicy); |
880 | > | if (myPolicy == "MIX") { |
881 | > | cp = MIX_CUTOFF_POLICY; |
882 | > | } else { |
883 | > | if (myPolicy == "MAX") { |
884 | > | cp = MAX_CUTOFF_POLICY; |
885 | > | } else { |
886 | > | if (myPolicy == "TRADITIONAL") { |
887 | > | cp = TRADITIONAL_CUTOFF_POLICY; |
888 | > | } else { |
889 | > | // throw error |
890 | > | sprintf( painCave.errMsg, |
891 | > | "SimInfo error: Unknown cutoffPolicy. (Input file specified %s .)\n\tcutoffPolicy must be one of: \"Mix\", \"Max\", or \"Traditional\".", myPolicy.c_str() ); |
892 | > | painCave.isFatal = 1; |
893 | > | simError(); |
894 | > | } |
895 | > | } |
896 | > | } |
897 | > | } |
898 | > | |
899 | > | |
900 | > | if (simParams_->haveSkinThickness()) { |
901 | > | double skinThickness = simParams_->getSkinThickness(); |
902 | > | } |
903 | > | |
904 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist, &cp); |
905 | > | // also send cutoff notification to electrostatics |
906 | > | setElectrostaticCutoffRadius(&rcut_, &rsw_); |
907 | > | } |
908 | > | |
909 | > | void SimInfo::setupElectrostaticSummationMethod( int isError ) { |
910 | > | |
911 | > | int errorOut; |
912 | > | int esm = NONE; |
913 | > | int sm = UNDAMPED; |
914 | > | double alphaVal; |
915 | > | double dielectric; |
916 | > | |
917 | > | errorOut = isError; |
918 | > | alphaVal = simParams_->getDampingAlpha(); |
919 | > | dielectric = simParams_->getDielectric(); |
920 | > | |
921 | > | if (simParams_->haveElectrostaticSummationMethod()) { |
922 | > | std::string myMethod = simParams_->getElectrostaticSummationMethod(); |
923 | > | toUpper(myMethod); |
924 | > | if (myMethod == "NONE") { |
925 | > | esm = NONE; |
926 | > | } else { |
927 | > | if (myMethod == "SWITCHING_FUNCTION") { |
928 | > | esm = SWITCHING_FUNCTION; |
929 | > | } else { |
930 | > | if (myMethod == "SHIFTED_POTENTIAL") { |
931 | > | esm = SHIFTED_POTENTIAL; |
932 | > | } else { |
933 | > | if (myMethod == "SHIFTED_FORCE") { |
934 | > | esm = SHIFTED_FORCE; |
935 | > | } else { |
936 | > | if (myMethod == "REACTION_FIELD") { |
937 | > | esm = REACTION_FIELD; |
938 | > | } else { |
939 | > | // throw error |
940 | > | sprintf( painCave.errMsg, |
941 | > | "SimInfo error: Unknown electrostaticSummationMethod. (Input file specified %s .)\n\telectrostaticSummationMethod must be one of: \"none\", \"shifted_potential\", \"shifted_force\", or \"reaction_field\".", myMethod.c_str() ); |
942 | > | painCave.isFatal = 1; |
943 | > | simError(); |
944 | > | } |
945 | > | } |
946 | > | } |
947 | > | } |
948 | > | } |
949 | > | } |
950 | > | |
951 | > | if (simParams_->haveElectrostaticScreeningMethod()) { |
952 | > | std::string myScreen = simParams_->getElectrostaticScreeningMethod(); |
953 | > | toUpper(myScreen); |
954 | > | if (myScreen == "UNDAMPED") { |
955 | > | sm = UNDAMPED; |
956 | > | } else { |
957 | > | if (myScreen == "DAMPED") { |
958 | > | sm = DAMPED; |
959 | > | if (!simParams_->haveDampingAlpha()) { |
960 | > | //throw error |
961 | > | sprintf( painCave.errMsg, |
962 | > | "SimInfo warning: dampingAlpha was not specified in the input file. A default value of %f (1/ang) will be used.", alphaVal); |
963 | > | painCave.isFatal = 0; |
964 | > | simError(); |
965 | > | } |
966 | > | } else { |
967 | > | // throw error |
968 | > | sprintf( painCave.errMsg, |
969 | > | "SimInfo error: Unknown electrostaticScreeningMethod. (Input file specified %s .)\n\telectrostaticScreeningMethod must be one of: \"undamped\" or \"damped\".", myScreen.c_str() ); |
970 | > | painCave.isFatal = 1; |
971 | > | simError(); |
972 | > | } |
973 | > | } |
974 | > | } |
975 | > | |
976 | > | // let's pass some summation method variables to fortran |
977 | > | setElectrostaticSummationMethod( &esm ); |
978 | > | setScreeningMethod( &sm ); |
979 | > | setDampingAlpha( &alphaVal ); |
980 | > | setReactionFieldDielectric( &dielectric ); |
981 | > | initFortranFF( &esm, &errorOut ); |
982 | > | } |
983 | > | |
984 | > | void SimInfo::setupSwitchingFunction() { |
985 | > | int ft = CUBIC; |
986 | > | |
987 | > | if (simParams_->haveSwitchingFunctionType()) { |
988 | > | std::string funcType = simParams_->getSwitchingFunctionType(); |
989 | > | toUpper(funcType); |
990 | > | if (funcType == "CUBIC") { |
991 | > | ft = CUBIC; |
992 | > | } else { |
993 | > | if (funcType == "FIFTH_ORDER_POLYNOMIAL") { |
994 | > | ft = FIFTH_ORDER_POLY; |
995 | > | } else { |
996 | > | // throw error |
997 | > | sprintf( painCave.errMsg, |
998 | > | "SimInfo error: Unknown switchingFunctionType. (Input file specified %s .)\n\tswitchingFunctionType must be one of: \"cubic\" or \"fifth_order_polynomial\".", funcType.c_str() ); |
999 | > | painCave.isFatal = 1; |
1000 | > | simError(); |
1001 | > | } |
1002 | > | } |
1003 | > | } |
1004 | > | |
1005 | > | // send switching function notification to switcheroo |
1006 | > | setFunctionType(&ft); |
1007 | > | |
1008 | > | } |
1009 | > | |
1010 | > | void SimInfo::addProperty(GenericData* genData) { |
1011 | > | properties_.addProperty(genData); |
1012 | > | } |
1013 | > | |
1014 | > | void SimInfo::removeProperty(const std::string& propName) { |
1015 | > | properties_.removeProperty(propName); |
1016 | > | } |
1017 | > | |
1018 | > | void SimInfo::clearProperties() { |
1019 | > | properties_.clearProperties(); |
1020 | > | } |
1021 | > | |
1022 | > | std::vector<std::string> SimInfo::getPropertyNames() { |
1023 | > | return properties_.getPropertyNames(); |
1024 | > | } |
1025 | > | |
1026 | > | std::vector<GenericData*> SimInfo::getProperties() { |
1027 | > | return properties_.getProperties(); |
1028 | > | } |
1029 | > | |
1030 | > | GenericData* SimInfo::getPropertyByName(const std::string& propName) { |
1031 | > | return properties_.getPropertyByName(propName); |
1032 | > | } |
1033 | > | |
1034 | > | void SimInfo::setSnapshotManager(SnapshotManager* sman) { |
1035 | > | if (sman_ == sman) { |
1036 | > | return; |
1037 | > | } |
1038 | > | delete sman_; |
1039 | > | sman_ = sman; |
1040 | > | |
1041 | > | Molecule* mol; |
1042 | > | RigidBody* rb; |
1043 | > | Atom* atom; |
1044 | > | SimInfo::MoleculeIterator mi; |
1045 | > | Molecule::RigidBodyIterator rbIter; |
1046 | > | Molecule::AtomIterator atomIter;; |
1047 | > | |
1048 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
1049 | > | |
1050 | > | for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { |
1051 | > | atom->setSnapshotManager(sman_); |
1052 | > | } |
1053 | > | |
1054 | > | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
1055 | > | rb->setSnapshotManager(sman_); |
1056 | > | } |
1057 | > | } |
1058 | > | |
1059 | > | } |
1060 | > | |
1061 | > | Vector3d SimInfo::getComVel(){ |
1062 | > | SimInfo::MoleculeIterator i; |
1063 | > | Molecule* mol; |
1064 | > | |
1065 | > | Vector3d comVel(0.0); |
1066 | > | double totalMass = 0.0; |
1067 | > | |
1068 | > | |
1069 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1070 | > | double mass = mol->getMass(); |
1071 | > | totalMass += mass; |
1072 | > | comVel += mass * mol->getComVel(); |
1073 | > | } |
1074 | > | |
1075 | > | #ifdef IS_MPI |
1076 | > | double tmpMass = totalMass; |
1077 | > | Vector3d tmpComVel(comVel); |
1078 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1079 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1080 | > | #endif |
1081 | > | |
1082 | > | comVel /= totalMass; |
1083 | > | |
1084 | > | return comVel; |
1085 | > | } |
1086 | > | |
1087 | > | Vector3d SimInfo::getCom(){ |
1088 | > | SimInfo::MoleculeIterator i; |
1089 | > | Molecule* mol; |
1090 | > | |
1091 | > | Vector3d com(0.0); |
1092 | > | double totalMass = 0.0; |
1093 | > | |
1094 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1095 | > | double mass = mol->getMass(); |
1096 | > | totalMass += mass; |
1097 | > | com += mass * mol->getCom(); |
1098 | > | } |
1099 | > | |
1100 | > | #ifdef IS_MPI |
1101 | > | double tmpMass = totalMass; |
1102 | > | Vector3d tmpCom(com); |
1103 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1104 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1105 | > | #endif |
1106 | > | |
1107 | > | com /= totalMass; |
1108 | > | |
1109 | > | return com; |
1110 | > | |
1111 | > | } |
1112 | > | |
1113 | > | std::ostream& operator <<(std::ostream& o, SimInfo& info) { |
1114 | > | |
1115 | > | return o; |
1116 | > | } |
1117 | > | |
1118 | > | |
1119 | > | /* |
1120 | > | Returns center of mass and center of mass velocity in one function call. |
1121 | > | */ |
1122 | > | |
1123 | > | void SimInfo::getComAll(Vector3d &com, Vector3d &comVel){ |
1124 | > | SimInfo::MoleculeIterator i; |
1125 | > | Molecule* mol; |
1126 | > | |
1127 | > | |
1128 | > | double totalMass = 0.0; |
1129 | > | |
1130 | > | |
1131 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1132 | > | double mass = mol->getMass(); |
1133 | > | totalMass += mass; |
1134 | > | com += mass * mol->getCom(); |
1135 | > | comVel += mass * mol->getComVel(); |
1136 | > | } |
1137 | > | |
1138 | > | #ifdef IS_MPI |
1139 | > | double tmpMass = totalMass; |
1140 | > | Vector3d tmpCom(com); |
1141 | > | Vector3d tmpComVel(comVel); |
1142 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1143 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1144 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1145 | > | #endif |
1146 | > | |
1147 | > | com /= totalMass; |
1148 | > | comVel /= totalMass; |
1149 | > | } |
1150 | > | |
1151 | > | /* |
1152 | > | Return intertia tensor for entire system and angular momentum Vector. |
1153 | > | |
1154 | > | |
1155 | > | [ Ixx -Ixy -Ixz ] |
1156 | > | J =| -Iyx Iyy -Iyz | |
1157 | > | [ -Izx -Iyz Izz ] |
1158 | > | */ |
1159 | > | |
1160 | > | void SimInfo::getInertiaTensor(Mat3x3d &inertiaTensor, Vector3d &angularMomentum){ |
1161 | > | |
1162 | > | |
1163 | > | double xx = 0.0; |
1164 | > | double yy = 0.0; |
1165 | > | double zz = 0.0; |
1166 | > | double xy = 0.0; |
1167 | > | double xz = 0.0; |
1168 | > | double yz = 0.0; |
1169 | > | Vector3d com(0.0); |
1170 | > | Vector3d comVel(0.0); |
1171 | > | |
1172 | > | getComAll(com, comVel); |
1173 | > | |
1174 | > | SimInfo::MoleculeIterator i; |
1175 | > | Molecule* mol; |
1176 | > | |
1177 | > | Vector3d thisq(0.0); |
1178 | > | Vector3d thisv(0.0); |
1179 | > | |
1180 | > | double thisMass = 0.0; |
1181 | > | |
1182 | > | |
1183 | > | |
1184 | > | |
1185 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1186 | > | |
1187 | > | thisq = mol->getCom()-com; |
1188 | > | thisv = mol->getComVel()-comVel; |
1189 | > | thisMass = mol->getMass(); |
1190 | > | // Compute moment of intertia coefficients. |
1191 | > | xx += thisq[0]*thisq[0]*thisMass; |
1192 | > | yy += thisq[1]*thisq[1]*thisMass; |
1193 | > | zz += thisq[2]*thisq[2]*thisMass; |
1194 | > | |
1195 | > | // compute products of intertia |
1196 | > | xy += thisq[0]*thisq[1]*thisMass; |
1197 | > | xz += thisq[0]*thisq[2]*thisMass; |
1198 | > | yz += thisq[1]*thisq[2]*thisMass; |
1199 | > | |
1200 | > | angularMomentum += cross( thisq, thisv ) * thisMass; |
1201 | > | |
1202 | > | } |
1203 | > | |
1204 | > | |
1205 | > | inertiaTensor(0,0) = yy + zz; |
1206 | > | inertiaTensor(0,1) = -xy; |
1207 | > | inertiaTensor(0,2) = -xz; |
1208 | > | inertiaTensor(1,0) = -xy; |
1209 | > | inertiaTensor(1,1) = xx + zz; |
1210 | > | inertiaTensor(1,2) = -yz; |
1211 | > | inertiaTensor(2,0) = -xz; |
1212 | > | inertiaTensor(2,1) = -yz; |
1213 | > | inertiaTensor(2,2) = xx + yy; |
1214 | > | |
1215 | > | #ifdef IS_MPI |
1216 | > | Mat3x3d tmpI(inertiaTensor); |
1217 | > | Vector3d tmpAngMom; |
1218 | > | MPI_Allreduce(tmpI.getArrayPointer(), inertiaTensor.getArrayPointer(),9,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1219 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1220 | > | #endif |
1221 | > | |
1222 | > | return; |
1223 | > | } |
1224 | > | |
1225 | > | //Returns the angular momentum of the system |
1226 | > | Vector3d SimInfo::getAngularMomentum(){ |
1227 | > | |
1228 | > | Vector3d com(0.0); |
1229 | > | Vector3d comVel(0.0); |
1230 | > | Vector3d angularMomentum(0.0); |
1231 | > | |
1232 | > | getComAll(com,comVel); |
1233 | > | |
1234 | > | SimInfo::MoleculeIterator i; |
1235 | > | Molecule* mol; |
1236 | > | |
1237 | > | Vector3d thisr(0.0); |
1238 | > | Vector3d thisp(0.0); |
1239 | > | |
1240 | > | double thisMass; |
1241 | > | |
1242 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
1243 | > | thisMass = mol->getMass(); |
1244 | > | thisr = mol->getCom()-com; |
1245 | > | thisp = (mol->getComVel()-comVel)*thisMass; |
1246 | > | |
1247 | > | angularMomentum += cross( thisr, thisp ); |
1248 | > | |
1249 | > | } |
1250 | > | |
1251 | > | #ifdef IS_MPI |
1252 | > | Vector3d tmpAngMom; |
1253 | > | MPI_Allreduce(tmpAngMom.getArrayPointer(), angularMomentum.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
1254 | > | #endif |
1255 | > | |
1256 | > | return angularMomentum; |
1257 | > | } |
1258 | > | |
1259 | > | |
1260 | > | }//end namespace oopse |
1261 | > |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |