# | 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 "SimInfo.hpp" |
53 | < | #define __C |
54 | < | #include "fSimulation.h" |
55 | < | #include "simError.h" |
52 | > | #include "brains/SimInfo.hpp" |
53 | > | #include "math/Vector3.hpp" |
54 | > | #include "primitives/Molecule.hpp" |
55 | > | #include "UseTheForce/doForces_interface.h" |
56 | > | #include "UseTheForce/notifyCutoffs_interface.h" |
57 | > | #include "utils/MemoryUtils.hpp" |
58 | > | #include "utils/simError.h" |
59 | > | #include "selection/SelectionManager.hpp" |
60 | ||
61 | < | #include "fortranWrappers.hpp" |
61 | > | #ifdef IS_MPI |
62 | > | #include "UseTheForce/mpiComponentPlan.h" |
63 | > | #include "UseTheForce/DarkSide/simParallel_interface.h" |
64 | > | #endif |
65 | ||
66 | < | #include "MatVec3.h" |
66 | > | namespace oopse { |
67 | ||
68 | < | #ifdef IS_MPI |
69 | < | #include "mpiSimulation.hpp" |
70 | < | #endif |
68 | > | SimInfo::SimInfo(std::vector<std::pair<MoleculeStamp*, int> >& molStampPairs, |
69 | > | ForceField* ff, Globals* simParams) : |
70 | > | forceField_(ff), simParams_(simParams), |
71 | > | ndf_(0), ndfRaw_(0), ndfTrans_(0), nZconstraint_(0), |
72 | > | nGlobalMols_(0), nGlobalAtoms_(0), nGlobalCutoffGroups_(0), |
73 | > | nGlobalIntegrableObjects_(0), nGlobalRigidBodies_(0), |
74 | > | nAtoms_(0), nBonds_(0), nBends_(0), nTorsions_(0), nRigidBodies_(0), |
75 | > | nIntegrableObjects_(0), nCutoffGroups_(0), nConstraints_(0), |
76 | > | sman_(NULL), fortranInitialized_(false) { |
77 | ||
78 | < | inline double roundMe( double x ){ |
79 | < | return ( x >= 0 ) ? floor( x + 0.5 ) : ceil( x - 0.5 ); |
80 | < | } |
81 | < | |
82 | < | inline double min( double a, double b ){ |
83 | < | return (a < b ) ? a : b; |
84 | < | } |
78 | > | |
79 | > | std::vector<std::pair<MoleculeStamp*, int> >::iterator i; |
80 | > | MoleculeStamp* molStamp; |
81 | > | int nMolWithSameStamp; |
82 | > | int nCutoffAtoms = 0; // number of atoms belong to cutoff groups |
83 | > | int nGroups = 0; //total cutoff groups defined in meta-data file |
84 | > | CutoffGroupStamp* cgStamp; |
85 | > | RigidBodyStamp* rbStamp; |
86 | > | int nRigidAtoms = 0; |
87 | > | |
88 | > | for (i = molStampPairs.begin(); i !=molStampPairs.end(); ++i) { |
89 | > | molStamp = i->first; |
90 | > | nMolWithSameStamp = i->second; |
91 | > | |
92 | > | addMoleculeStamp(molStamp, nMolWithSameStamp); |
93 | ||
94 | < | SimInfo* currentInfo; |
94 | > | //calculate atoms in molecules |
95 | > | nGlobalAtoms_ += molStamp->getNAtoms() *nMolWithSameStamp; |
96 | ||
31 | – | SimInfo::SimInfo(){ |
97 | ||
98 | < | n_constraints = 0; |
99 | < | nZconstraints = 0; |
100 | < | n_oriented = 0; |
101 | < | n_dipoles = 0; |
102 | < | ndf = 0; |
103 | < | ndfRaw = 0; |
104 | < | nZconstraints = 0; |
105 | < | the_integrator = NULL; |
41 | < | setTemp = 0; |
42 | < | thermalTime = 0.0; |
43 | < | currentTime = 0.0; |
44 | < | rCut = 0.0; |
45 | < | rSw = 0.0; |
98 | > | //calculate atoms in cutoff groups |
99 | > | int nAtomsInGroups = 0; |
100 | > | int nCutoffGroupsInStamp = molStamp->getNCutoffGroups(); |
101 | > | |
102 | > | for (int j=0; j < nCutoffGroupsInStamp; j++) { |
103 | > | cgStamp = molStamp->getCutoffGroup(j); |
104 | > | nAtomsInGroups += cgStamp->getNMembers(); |
105 | > | } |
106 | ||
107 | < | haveRcut = 0; |
108 | < | haveRsw = 0; |
49 | < | boxIsInit = 0; |
50 | < | |
51 | < | resetTime = 1e99; |
107 | > | nGroups += nCutoffGroupsInStamp * nMolWithSameStamp; |
108 | > | nCutoffAtoms += nAtomsInGroups * nMolWithSameStamp; |
109 | ||
110 | < | orthoRhombic = 0; |
111 | < | orthoTolerance = 1E-6; |
112 | < | useInitXSstate = true; |
110 | > | //calculate atoms in rigid bodies |
111 | > | int nAtomsInRigidBodies = 0; |
112 | > | int nRigidBodiesInStamp = molStamp->getNRigidBodies(); |
113 | > | |
114 | > | for (int j=0; j < nRigidBodiesInStamp; j++) { |
115 | > | rbStamp = molStamp->getRigidBody(j); |
116 | > | nAtomsInRigidBodies += rbStamp->getNMembers(); |
117 | > | } |
118 | ||
119 | < | usePBC = 0; |
120 | < | useLJ = 0; |
121 | < | useSticky = 0; |
122 | < | useCharges = 0; |
61 | < | useDipoles = 0; |
62 | < | useReactionField = 0; |
63 | < | useGB = 0; |
64 | < | useEAM = 0; |
65 | < | useSolidThermInt = 0; |
66 | < | useLiquidThermInt = 0; |
119 | > | nGlobalRigidBodies_ += nRigidBodiesInStamp * nMolWithSameStamp; |
120 | > | nRigidAtoms += nAtomsInRigidBodies * nMolWithSameStamp; |
121 | > | |
122 | > | } |
123 | ||
124 | < | haveCutoffGroups = false; |
124 | > | //every free atom (atom does not belong to cutoff groups) is a cutoff group |
125 | > | //therefore the total number of cutoff groups in the system is equal to |
126 | > | //the total number of atoms minus number of atoms belong to cutoff group defined in meta-data |
127 | > | //file plus the number of cutoff groups defined in meta-data file |
128 | > | nGlobalCutoffGroups_ = nGlobalAtoms_ - nCutoffAtoms + nGroups; |
129 | ||
130 | < | excludes = Exclude::Instance(); |
130 | > | //every free atom (atom does not belong to rigid bodies) is an integrable object |
131 | > | //therefore the total number of integrable objects in the system is equal to |
132 | > | //the total number of atoms minus number of atoms belong to rigid body defined in meta-data |
133 | > | //file plus the number of rigid bodies defined in meta-data file |
134 | > | nGlobalIntegrableObjects_ = nGlobalAtoms_ - nRigidAtoms + nGlobalRigidBodies_; |
135 | ||
136 | < | myConfiguration = new SimState(); |
136 | > | nGlobalMols_ = molStampIds_.size(); |
137 | ||
138 | < | has_minimizer = false; |
139 | < | the_minimizer =NULL; |
138 | > | #ifdef IS_MPI |
139 | > | molToProcMap_.resize(nGlobalMols_); |
140 | > | #endif |
141 | ||
142 | < | ngroup = 0; |
142 | > | } |
143 | ||
144 | < | wrapMeSimInfo( this ); |
144 | > | SimInfo::~SimInfo() { |
145 | > | std::map<int, Molecule*>::iterator i; |
146 | > | for (i = molecules_.begin(); i != molecules_.end(); ++i) { |
147 | > | delete i->second; |
148 | > | } |
149 | > | molecules_.clear(); |
150 | > | |
151 | > | MemoryUtils::deletePointers(moleculeStamps_); |
152 | > | |
153 | > | delete sman_; |
154 | > | delete simParams_; |
155 | > | delete forceField_; |
156 | } | |
157 | ||
158 | + | int SimInfo::getNGlobalConstraints() { |
159 | + | int nGlobalConstraints; |
160 | + | #ifdef IS_MPI |
161 | + | MPI_Allreduce(&nConstraints_, &nGlobalConstraints, 1, MPI_INT, MPI_SUM, |
162 | + | MPI_COMM_WORLD); |
163 | + | #else |
164 | + | nGlobalConstraints = nConstraints_; |
165 | + | #endif |
166 | + | return nGlobalConstraints; |
167 | + | } |
168 | ||
169 | < | SimInfo::~SimInfo(){ |
169 | > | bool SimInfo::addMolecule(Molecule* mol) { |
170 | > | MoleculeIterator i; |
171 | ||
172 | < | delete myConfiguration; |
172 | > | i = molecules_.find(mol->getGlobalIndex()); |
173 | > | if (i == molecules_.end() ) { |
174 | ||
175 | < | map<string, GenericData*>::iterator i; |
176 | < | |
177 | < | for(i = properties.begin(); i != properties.end(); i++) |
178 | < | delete (*i).second; |
175 | > | molecules_.insert(std::make_pair(mol->getGlobalIndex(), mol)); |
176 | > | |
177 | > | nAtoms_ += mol->getNAtoms(); |
178 | > | nBonds_ += mol->getNBonds(); |
179 | > | nBends_ += mol->getNBends(); |
180 | > | nTorsions_ += mol->getNTorsions(); |
181 | > | nRigidBodies_ += mol->getNRigidBodies(); |
182 | > | nIntegrableObjects_ += mol->getNIntegrableObjects(); |
183 | > | nCutoffGroups_ += mol->getNCutoffGroups(); |
184 | > | nConstraints_ += mol->getNConstraintPairs(); |
185 | ||
186 | + | addExcludePairs(mol); |
187 | + | |
188 | + | return true; |
189 | + | } else { |
190 | + | return false; |
191 | + | } |
192 | } | |
193 | ||
194 | < | void SimInfo::setBox(double newBox[3]) { |
195 | < | |
196 | < | int i, j; |
97 | < | double tempMat[3][3]; |
194 | > | bool SimInfo::removeMolecule(Molecule* mol) { |
195 | > | MoleculeIterator i; |
196 | > | i = molecules_.find(mol->getGlobalIndex()); |
197 | ||
198 | < | for(i=0; i<3; i++) |
100 | < | for (j=0; j<3; j++) tempMat[i][j] = 0.0;; |
198 | > | if (i != molecules_.end() ) { |
199 | ||
200 | < | tempMat[0][0] = newBox[0]; |
201 | < | tempMat[1][1] = newBox[1]; |
202 | < | tempMat[2][2] = newBox[2]; |
200 | > | assert(mol == i->second); |
201 | > | |
202 | > | nAtoms_ -= mol->getNAtoms(); |
203 | > | nBonds_ -= mol->getNBonds(); |
204 | > | nBends_ -= mol->getNBends(); |
205 | > | nTorsions_ -= mol->getNTorsions(); |
206 | > | nRigidBodies_ -= mol->getNRigidBodies(); |
207 | > | nIntegrableObjects_ -= mol->getNIntegrableObjects(); |
208 | > | nCutoffGroups_ -= mol->getNCutoffGroups(); |
209 | > | nConstraints_ -= mol->getNConstraintPairs(); |
210 | ||
211 | < | setBoxM( tempMat ); |
211 | > | removeExcludePairs(mol); |
212 | > | molecules_.erase(mol->getGlobalIndex()); |
213 | ||
214 | < | } |
214 | > | delete mol; |
215 | > | |
216 | > | return true; |
217 | > | } else { |
218 | > | return false; |
219 | > | } |
220 | ||
110 | – | void SimInfo::setBoxM( double theBox[3][3] ){ |
111 | – | |
112 | – | int i, j; |
113 | – | double FortranHmat[9]; // to preserve compatibility with Fortran the |
114 | – | // ordering in the array is as follows: |
115 | – | // [ 0 3 6 ] |
116 | – | // [ 1 4 7 ] |
117 | – | // [ 2 5 8 ] |
118 | – | double FortranHmatInv[9]; // the inverted Hmat (for Fortran); |
221 | ||
222 | < | if( !boxIsInit ) boxIsInit = 1; |
222 | > | } |
223 | ||
224 | < | for(i=0; i < 3; i++) |
225 | < | for (j=0; j < 3; j++) Hmat[i][j] = theBox[i][j]; |
226 | < | |
227 | < | calcBoxL(); |
228 | < | calcHmatInv(); |
224 | > | |
225 | > | Molecule* SimInfo::beginMolecule(MoleculeIterator& i) { |
226 | > | i = molecules_.begin(); |
227 | > | return i == molecules_.end() ? NULL : i->second; |
228 | > | } |
229 | ||
230 | < | for(i=0; i < 3; i++) { |
231 | < | for (j=0; j < 3; j++) { |
232 | < | FortranHmat[3*j + i] = Hmat[i][j]; |
131 | < | FortranHmatInv[3*j + i] = HmatInv[i][j]; |
132 | < | } |
133 | < | } |
134 | < | |
135 | < | setFortranBoxSize(FortranHmat, FortranHmatInv, &orthoRhombic); |
136 | < | |
230 | > | Molecule* SimInfo::nextMolecule(MoleculeIterator& i) { |
231 | > | ++i; |
232 | > | return i == molecules_.end() ? NULL : i->second; |
233 | } | |
138 | – | |
234 | ||
140 | – | void SimInfo::getBoxM (double theBox[3][3]) { |
235 | ||
236 | < | int i, j; |
237 | < | for(i=0; i<3; i++) |
238 | < | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]; |
239 | < | } |
236 | > | void SimInfo::calcNdf() { |
237 | > | int ndf_local; |
238 | > | MoleculeIterator i; |
239 | > | std::vector<StuntDouble*>::iterator j; |
240 | > | Molecule* mol; |
241 | > | StuntDouble* integrableObject; |
242 | ||
243 | + | ndf_local = 0; |
244 | + | |
245 | + | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
246 | + | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
247 | + | integrableObject = mol->nextIntegrableObject(j)) { |
248 | ||
249 | < | void SimInfo::scaleBox(double scale) { |
149 | < | double theBox[3][3]; |
150 | < | int i, j; |
249 | > | ndf_local += 3; |
250 | ||
251 | < | // cerr << "Scaling box by " << scale << "\n"; |
251 | > | if (integrableObject->isDirectional()) { |
252 | > | if (integrableObject->isLinear()) { |
253 | > | ndf_local += 2; |
254 | > | } else { |
255 | > | ndf_local += 3; |
256 | > | } |
257 | > | } |
258 | > | |
259 | > | }//end for (integrableObject) |
260 | > | }// end for (mol) |
261 | > | |
262 | > | // n_constraints is local, so subtract them on each processor |
263 | > | ndf_local -= nConstraints_; |
264 | ||
265 | < | for(i=0; i<3; i++) |
266 | < | for (j=0; j<3; j++) theBox[i][j] = Hmat[i][j]*scale; |
265 | > | #ifdef IS_MPI |
266 | > | MPI_Allreduce(&ndf_local,&ndf_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
267 | > | #else |
268 | > | ndf_ = ndf_local; |
269 | > | #endif |
270 | ||
271 | < | setBoxM(theBox); |
271 | > | // nZconstraints_ is global, as are the 3 COM translations for the |
272 | > | // entire system: |
273 | > | ndf_ = ndf_ - 3 - nZconstraint_; |
274 | ||
275 | } | |
276 | ||
277 | < | void SimInfo::calcHmatInv( void ) { |
278 | < | |
163 | < | int oldOrtho; |
164 | < | int i,j; |
165 | < | double smallDiag; |
166 | < | double tol; |
167 | < | double sanity[3][3]; |
277 | > | void SimInfo::calcNdfRaw() { |
278 | > | int ndfRaw_local; |
279 | ||
280 | < | invertMat3( Hmat, HmatInv ); |
280 | > | MoleculeIterator i; |
281 | > | std::vector<StuntDouble*>::iterator j; |
282 | > | Molecule* mol; |
283 | > | StuntDouble* integrableObject; |
284 | ||
285 | < | // check to see if Hmat is orthorhombic |
286 | < | |
287 | < | oldOrtho = orthoRhombic; |
285 | > | // Raw degrees of freedom that we have to set |
286 | > | ndfRaw_local = 0; |
287 | > | |
288 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
289 | > | for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL; |
290 | > | integrableObject = mol->nextIntegrableObject(j)) { |
291 | ||
292 | < | smallDiag = fabs(Hmat[0][0]); |
176 | < | if(smallDiag > fabs(Hmat[1][1])) smallDiag = fabs(Hmat[1][1]); |
177 | < | if(smallDiag > fabs(Hmat[2][2])) smallDiag = fabs(Hmat[2][2]); |
178 | < | tol = smallDiag * orthoTolerance; |
292 | > | ndfRaw_local += 3; |
293 | ||
294 | < | orthoRhombic = 1; |
295 | < | |
296 | < | for (i = 0; i < 3; i++ ) { |
297 | < | for (j = 0 ; j < 3; j++) { |
298 | < | if (i != j) { |
299 | < | if (orthoRhombic) { |
300 | < | if ( fabs(Hmat[i][j]) >= tol) orthoRhombic = 0; |
301 | < | } |
302 | < | } |
294 | > | if (integrableObject->isDirectional()) { |
295 | > | if (integrableObject->isLinear()) { |
296 | > | ndfRaw_local += 2; |
297 | > | } else { |
298 | > | ndfRaw_local += 3; |
299 | > | } |
300 | > | } |
301 | > | |
302 | > | } |
303 | } | |
190 | – | } |
191 | – | |
192 | – | if( oldOrtho != orthoRhombic ){ |
304 | ||
305 | < | if( orthoRhombic ) { |
306 | < | sprintf( painCave.errMsg, |
307 | < | "OOPSE is switching from the default Non-Orthorhombic\n" |
308 | < | "\tto the faster Orthorhombic periodic boundary computations.\n" |
309 | < | "\tThis is usually a good thing, but if you wan't the\n" |
199 | < | "\tNon-Orthorhombic computations, make the orthoBoxTolerance\n" |
200 | < | "\tvariable ( currently set to %G ) smaller.\n", |
201 | < | orthoTolerance); |
202 | < | painCave.severity = OOPSE_INFO; |
203 | < | simError(); |
204 | < | } |
205 | < | else { |
206 | < | sprintf( painCave.errMsg, |
207 | < | "OOPSE is switching from the faster Orthorhombic to the more\n" |
208 | < | "\tflexible Non-Orthorhombic periodic boundary computations.\n" |
209 | < | "\tThis is usually because the box has deformed under\n" |
210 | < | "\tNPTf integration. If you wan't to live on the edge with\n" |
211 | < | "\tthe Orthorhombic computations, make the orthoBoxTolerance\n" |
212 | < | "\tvariable ( currently set to %G ) larger.\n", |
213 | < | orthoTolerance); |
214 | < | painCave.severity = OOPSE_WARNING; |
215 | < | simError(); |
216 | < | } |
217 | < | } |
305 | > | #ifdef IS_MPI |
306 | > | MPI_Allreduce(&ndfRaw_local,&ndfRaw_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
307 | > | #else |
308 | > | ndfRaw_ = ndfRaw_local; |
309 | > | #endif |
310 | } | |
311 | ||
312 | < | void SimInfo::calcBoxL( void ){ |
312 | > | void SimInfo::calcNdfTrans() { |
313 | > | int ndfTrans_local; |
314 | ||
315 | < | double dx, dy, dz, dsq; |
315 | > | ndfTrans_local = 3 * nIntegrableObjects_ - nConstraints_; |
316 | ||
224 | – | // boxVol = Determinant of Hmat |
317 | ||
318 | < | boxVol = matDet3( Hmat ); |
318 | > | #ifdef IS_MPI |
319 | > | MPI_Allreduce(&ndfTrans_local,&ndfTrans_,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
320 | > | #else |
321 | > | ndfTrans_ = ndfTrans_local; |
322 | > | #endif |
323 | ||
324 | < | // boxLx |
325 | < | |
326 | < | dx = Hmat[0][0]; dy = Hmat[1][0]; dz = Hmat[2][0]; |
231 | < | dsq = dx*dx + dy*dy + dz*dz; |
232 | < | boxL[0] = sqrt( dsq ); |
233 | < | //maxCutoff = 0.5 * boxL[0]; |
324 | > | ndfTrans_ = ndfTrans_ - 3 - nZconstraint_; |
325 | > | |
326 | > | } |
327 | ||
328 | < | // boxLy |
329 | < | |
330 | < | dx = Hmat[0][1]; dy = Hmat[1][1]; dz = Hmat[2][1]; |
331 | < | dsq = dx*dx + dy*dy + dz*dz; |
332 | < | boxL[1] = sqrt( dsq ); |
333 | < | //if( (0.5 * boxL[1]) < maxCutoff ) maxCutoff = 0.5 * boxL[1]; |
328 | > | void SimInfo::addExcludePairs(Molecule* mol) { |
329 | > | std::vector<Bond*>::iterator bondIter; |
330 | > | std::vector<Bend*>::iterator bendIter; |
331 | > | std::vector<Torsion*>::iterator torsionIter; |
332 | > | Bond* bond; |
333 | > | Bend* bend; |
334 | > | Torsion* torsion; |
335 | > | int a; |
336 | > | int b; |
337 | > | int c; |
338 | > | int d; |
339 | > | |
340 | > | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
341 | > | a = bond->getAtomA()->getGlobalIndex(); |
342 | > | b = bond->getAtomB()->getGlobalIndex(); |
343 | > | exclude_.addPair(a, b); |
344 | > | } |
345 | ||
346 | + | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
347 | + | a = bend->getAtomA()->getGlobalIndex(); |
348 | + | b = bend->getAtomB()->getGlobalIndex(); |
349 | + | c = bend->getAtomC()->getGlobalIndex(); |
350 | ||
351 | < | // boxLz |
352 | < | |
353 | < | dx = Hmat[0][2]; dy = Hmat[1][2]; dz = Hmat[2][2]; |
354 | < | dsq = dx*dx + dy*dy + dz*dz; |
247 | < | boxL[2] = sqrt( dsq ); |
248 | < | //if( (0.5 * boxL[2]) < maxCutoff ) maxCutoff = 0.5 * boxL[2]; |
351 | > | exclude_.addPair(a, b); |
352 | > | exclude_.addPair(a, c); |
353 | > | exclude_.addPair(b, c); |
354 | > | } |
355 | ||
356 | < | //calculate the max cutoff |
357 | < | maxCutoff = calcMaxCutOff(); |
358 | < | |
359 | < | checkCutOffs(); |
356 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
357 | > | a = torsion->getAtomA()->getGlobalIndex(); |
358 | > | b = torsion->getAtomB()->getGlobalIndex(); |
359 | > | c = torsion->getAtomC()->getGlobalIndex(); |
360 | > | d = torsion->getAtomD()->getGlobalIndex(); |
361 | ||
362 | + | exclude_.addPair(a, b); |
363 | + | exclude_.addPair(a, c); |
364 | + | exclude_.addPair(a, d); |
365 | + | exclude_.addPair(b, c); |
366 | + | exclude_.addPair(b, d); |
367 | + | exclude_.addPair(c, d); |
368 | + | } |
369 | + | |
370 | + | |
371 | } | |
372 | ||
373 | + | void SimInfo::removeExcludePairs(Molecule* mol) { |
374 | + | std::vector<Bond*>::iterator bondIter; |
375 | + | std::vector<Bend*>::iterator bendIter; |
376 | + | std::vector<Torsion*>::iterator torsionIter; |
377 | + | Bond* bond; |
378 | + | Bend* bend; |
379 | + | Torsion* torsion; |
380 | + | int a; |
381 | + | int b; |
382 | + | int c; |
383 | + | int d; |
384 | + | |
385 | + | for (bond= mol->beginBond(bondIter); bond != NULL; bond = mol->nextBond(bondIter)) { |
386 | + | a = bond->getAtomA()->getGlobalIndex(); |
387 | + | b = bond->getAtomB()->getGlobalIndex(); |
388 | + | exclude_.removePair(a, b); |
389 | + | } |
390 | ||
391 | < | double SimInfo::calcMaxCutOff(){ |
391 | > | for (bend= mol->beginBend(bendIter); bend != NULL; bend = mol->nextBend(bendIter)) { |
392 | > | a = bend->getAtomA()->getGlobalIndex(); |
393 | > | b = bend->getAtomB()->getGlobalIndex(); |
394 | > | c = bend->getAtomC()->getGlobalIndex(); |
395 | ||
396 | < | double ri[3], rj[3], rk[3]; |
397 | < | double rij[3], rjk[3], rki[3]; |
398 | < | double minDist; |
396 | > | exclude_.removePair(a, b); |
397 | > | exclude_.removePair(a, c); |
398 | > | exclude_.removePair(b, c); |
399 | > | } |
400 | ||
401 | < | ri[0] = Hmat[0][0]; |
402 | < | ri[1] = Hmat[1][0]; |
403 | < | ri[2] = Hmat[2][0]; |
401 | > | for (torsion= mol->beginTorsion(torsionIter); torsion != NULL; torsion = mol->nextTorsion(torsionIter)) { |
402 | > | a = torsion->getAtomA()->getGlobalIndex(); |
403 | > | b = torsion->getAtomB()->getGlobalIndex(); |
404 | > | c = torsion->getAtomC()->getGlobalIndex(); |
405 | > | d = torsion->getAtomD()->getGlobalIndex(); |
406 | ||
407 | < | rj[0] = Hmat[0][1]; |
408 | < | rj[1] = Hmat[1][1]; |
409 | < | rj[2] = Hmat[2][1]; |
407 | > | exclude_.removePair(a, b); |
408 | > | exclude_.removePair(a, c); |
409 | > | exclude_.removePair(a, d); |
410 | > | exclude_.removePair(b, c); |
411 | > | exclude_.removePair(b, d); |
412 | > | exclude_.removePair(c, d); |
413 | > | } |
414 | ||
415 | < | rk[0] = Hmat[0][2]; |
273 | < | rk[1] = Hmat[1][2]; |
274 | < | rk[2] = Hmat[2][2]; |
275 | < | |
276 | < | crossProduct3(ri, rj, rij); |
277 | < | distXY = dotProduct3(rk,rij) / norm3(rij); |
415 | > | } |
416 | ||
279 | – | crossProduct3(rj,rk, rjk); |
280 | – | distYZ = dotProduct3(ri,rjk) / norm3(rjk); |
417 | ||
418 | < | crossProduct3(rk,ri, rki); |
419 | < | distZX = dotProduct3(rj,rki) / norm3(rki); |
418 | > | void SimInfo::addMoleculeStamp(MoleculeStamp* molStamp, int nmol) { |
419 | > | int curStampId; |
420 | ||
421 | < | minDist = min(min(distXY, distYZ), distZX); |
422 | < | return minDist/2; |
423 | < | |
421 | > | //index from 0 |
422 | > | curStampId = moleculeStamps_.size(); |
423 | > | |
424 | > | moleculeStamps_.push_back(molStamp); |
425 | > | molStampIds_.insert(molStampIds_.end(), nmol, curStampId); |
426 | } | |
427 | ||
428 | < | void SimInfo::wrapVector( double thePos[3] ){ |
428 | > | void SimInfo::update() { |
429 | ||
430 | < | int i; |
293 | < | double scaled[3]; |
430 | > | setupSimType(); |
431 | ||
432 | < | if( !orthoRhombic ){ |
433 | < | // calc the scaled coordinates. |
434 | < | |
432 | > | #ifdef IS_MPI |
433 | > | setupFortranParallel(); |
434 | > | #endif |
435 | ||
436 | < | matVecMul3(HmatInv, thePos, scaled); |
300 | < | |
301 | < | for(i=0; i<3; i++) |
302 | < | scaled[i] -= roundMe(scaled[i]); |
303 | < | |
304 | < | // calc the wrapped real coordinates from the wrapped scaled coordinates |
305 | < | |
306 | < | matVecMul3(Hmat, scaled, thePos); |
436 | > | setupFortranSim(); |
437 | ||
438 | < | } |
439 | < | else{ |
440 | < | // calc the scaled coordinates. |
438 | > | //setup fortran force field |
439 | > | /** @deprecate */ |
440 | > | int isError = 0; |
441 | > | initFortranFF( &fInfo_.SIM_uses_RF , &isError ); |
442 | > | if(isError){ |
443 | > | sprintf( painCave.errMsg, |
444 | > | "ForceField error: There was an error initializing the forceField in fortran.\n" ); |
445 | > | painCave.isFatal = 1; |
446 | > | simError(); |
447 | > | } |
448 | > | |
449 | ||
450 | < | for(i=0; i<3; i++) |
451 | < | scaled[i] = thePos[i]*HmatInv[i][i]; |
452 | < | |
453 | < | // wrap the scaled coordinates |
454 | < | |
455 | < | for(i=0; i<3; i++) |
456 | < | scaled[i] -= roundMe(scaled[i]); |
319 | < | |
320 | < | // calc the wrapped real coordinates from the wrapped scaled coordinates |
321 | < | |
322 | < | for(i=0; i<3; i++) |
323 | < | thePos[i] = scaled[i]*Hmat[i][i]; |
324 | < | } |
325 | < | |
450 | > | setupCutoff(); |
451 | > | |
452 | > | calcNdf(); |
453 | > | calcNdfRaw(); |
454 | > | calcNdfTrans(); |
455 | > | |
456 | > | fortranInitialized_ = true; |
457 | } | |
458 | ||
459 | + | std::set<AtomType*> SimInfo::getUniqueAtomTypes() { |
460 | + | SimInfo::MoleculeIterator mi; |
461 | + | Molecule* mol; |
462 | + | Molecule::AtomIterator ai; |
463 | + | Atom* atom; |
464 | + | std::set<AtomType*> atomTypes; |
465 | ||
466 | < | int SimInfo::getNDF(){ |
330 | < | int ndf_local; |
466 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
467 | ||
468 | < | ndf_local = 0; |
469 | < | |
470 | < | for(int i = 0; i < integrableObjects.size(); i++){ |
471 | < | ndf_local += 3; |
336 | < | if (integrableObjects[i]->isDirectional()) { |
337 | < | if (integrableObjects[i]->isLinear()) |
338 | < | ndf_local += 2; |
339 | < | else |
340 | < | ndf_local += 3; |
468 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
469 | > | atomTypes.insert(atom->getAtomType()); |
470 | > | } |
471 | > | |
472 | } | |
342 | – | } |
473 | ||
474 | < | // n_constraints is local, so subtract them on each processor: |
474 | > | return atomTypes; |
475 | > | } |
476 | ||
477 | < | ndf_local -= n_constraints; |
477 | > | void SimInfo::setupSimType() { |
478 | > | std::set<AtomType*>::iterator i; |
479 | > | std::set<AtomType*> atomTypes; |
480 | > | atomTypes = getUniqueAtomTypes(); |
481 | > | |
482 | > | int useLennardJones = 0; |
483 | > | int useElectrostatic = 0; |
484 | > | int useEAM = 0; |
485 | > | int useCharge = 0; |
486 | > | int useDirectional = 0; |
487 | > | int useDipole = 0; |
488 | > | int useGayBerne = 0; |
489 | > | int useSticky = 0; |
490 | > | int useShape = 0; |
491 | > | int useFLARB = 0; //it is not in AtomType yet |
492 | > | int useDirectionalAtom = 0; |
493 | > | int useElectrostatics = 0; |
494 | > | //usePBC and useRF are from simParams |
495 | > | int usePBC = simParams_->getPBC(); |
496 | > | int useRF = simParams_->getUseRF(); |
497 | ||
498 | < | #ifdef IS_MPI |
499 | < | MPI_Allreduce(&ndf_local,&ndf,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
500 | < | #else |
501 | < | ndf = ndf_local; |
498 | > | //loop over all of the atom types |
499 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
500 | > | useLennardJones |= (*i)->isLennardJones(); |
501 | > | useElectrostatic |= (*i)->isElectrostatic(); |
502 | > | useEAM |= (*i)->isEAM(); |
503 | > | useCharge |= (*i)->isCharge(); |
504 | > | useDirectional |= (*i)->isDirectional(); |
505 | > | useDipole |= (*i)->isDipole(); |
506 | > | useGayBerne |= (*i)->isGayBerne(); |
507 | > | useSticky |= (*i)->isSticky(); |
508 | > | useShape |= (*i)->isShape(); |
509 | > | } |
510 | > | |
511 | > | if (useSticky || useDipole || useGayBerne || useShape) { |
512 | > | useDirectionalAtom = 1; |
513 | > | } |
514 | > | |
515 | > | if (useCharge || useDipole) { |
516 | > | useElectrostatics = 1; |
517 | > | } |
518 | > | |
519 | > | #ifdef IS_MPI |
520 | > | int temp; |
521 | > | |
522 | > | temp = usePBC; |
523 | > | MPI_Allreduce(&temp, &usePBC, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
524 | > | |
525 | > | temp = useDirectionalAtom; |
526 | > | MPI_Allreduce(&temp, &useDirectionalAtom, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
527 | > | |
528 | > | temp = useLennardJones; |
529 | > | MPI_Allreduce(&temp, &useLennardJones, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
530 | > | |
531 | > | temp = useElectrostatics; |
532 | > | MPI_Allreduce(&temp, &useElectrostatics, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
533 | > | |
534 | > | temp = useCharge; |
535 | > | MPI_Allreduce(&temp, &useCharge, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
536 | > | |
537 | > | temp = useDipole; |
538 | > | MPI_Allreduce(&temp, &useDipole, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
539 | > | |
540 | > | temp = useSticky; |
541 | > | MPI_Allreduce(&temp, &useSticky, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
542 | > | |
543 | > | temp = useGayBerne; |
544 | > | MPI_Allreduce(&temp, &useGayBerne, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
545 | > | |
546 | > | temp = useEAM; |
547 | > | MPI_Allreduce(&temp, &useEAM, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
548 | > | |
549 | > | temp = useShape; |
550 | > | MPI_Allreduce(&temp, &useShape, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
551 | > | |
552 | > | temp = useFLARB; |
553 | > | MPI_Allreduce(&temp, &useFLARB, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
554 | > | |
555 | > | temp = useRF; |
556 | > | MPI_Allreduce(&temp, &useRF, 1, MPI_INT, MPI_LOR, MPI_COMM_WORLD); |
557 | > | |
558 | #endif | |
559 | ||
560 | < | // nZconstraints is global, as are the 3 COM translations for the |
561 | < | // entire system: |
560 | > | fInfo_.SIM_uses_PBC = usePBC; |
561 | > | fInfo_.SIM_uses_DirectionalAtoms = useDirectionalAtom; |
562 | > | fInfo_.SIM_uses_LennardJones = useLennardJones; |
563 | > | fInfo_.SIM_uses_Electrostatics = useElectrostatics; |
564 | > | fInfo_.SIM_uses_Charges = useCharge; |
565 | > | fInfo_.SIM_uses_Dipoles = useDipole; |
566 | > | fInfo_.SIM_uses_Sticky = useSticky; |
567 | > | fInfo_.SIM_uses_GayBerne = useGayBerne; |
568 | > | fInfo_.SIM_uses_EAM = useEAM; |
569 | > | fInfo_.SIM_uses_Shapes = useShape; |
570 | > | fInfo_.SIM_uses_FLARB = useFLARB; |
571 | > | fInfo_.SIM_uses_RF = useRF; |
572 | ||
573 | < | ndf = ndf - 3 - nZconstraints; |
573 | > | if( fInfo_.SIM_uses_Dipoles && fInfo_.SIM_uses_RF) { |
574 | ||
575 | < | return ndf; |
575 | > | if (simParams_->haveDielectric()) { |
576 | > | fInfo_.dielect = simParams_->getDielectric(); |
577 | > | } else { |
578 | > | sprintf(painCave.errMsg, |
579 | > | "SimSetup Error: No Dielectric constant was set.\n" |
580 | > | "\tYou are trying to use Reaction Field without" |
581 | > | "\tsetting a dielectric constant!\n"); |
582 | > | painCave.isFatal = 1; |
583 | > | simError(); |
584 | > | } |
585 | > | |
586 | > | } else { |
587 | > | fInfo_.dielect = 0.0; |
588 | > | } |
589 | > | |
590 | } | |
591 | ||
592 | < | int SimInfo::getNDFraw() { |
593 | < | int ndfRaw_local; |
592 | > | void SimInfo::setupFortranSim() { |
593 | > | int isError; |
594 | > | int nExclude; |
595 | > | std::vector<int> fortranGlobalGroupMembership; |
596 | > | |
597 | > | nExclude = exclude_.getSize(); |
598 | > | isError = 0; |
599 | ||
600 | < | // Raw degrees of freedom that we have to set |
601 | < | ndfRaw_local = 0; |
600 | > | //globalGroupMembership_ is filled by SimCreator |
601 | > | for (int i = 0; i < nGlobalAtoms_; i++) { |
602 | > | fortranGlobalGroupMembership.push_back(globalGroupMembership_[i] + 1); |
603 | > | } |
604 | ||
605 | < | for(int i = 0; i < integrableObjects.size(); i++){ |
606 | < | ndfRaw_local += 3; |
607 | < | if (integrableObjects[i]->isDirectional()) { |
608 | < | if (integrableObjects[i]->isLinear()) |
609 | < | ndfRaw_local += 2; |
610 | < | else |
611 | < | ndfRaw_local += 3; |
605 | > | //calculate mass ratio of cutoff group |
606 | > | std::vector<double> mfact; |
607 | > | SimInfo::MoleculeIterator mi; |
608 | > | Molecule* mol; |
609 | > | Molecule::CutoffGroupIterator ci; |
610 | > | CutoffGroup* cg; |
611 | > | Molecule::AtomIterator ai; |
612 | > | Atom* atom; |
613 | > | double totalMass; |
614 | > | |
615 | > | //to avoid memory reallocation, reserve enough space for mfact |
616 | > | mfact.reserve(getNCutoffGroups()); |
617 | > | |
618 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
619 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
620 | > | |
621 | > | totalMass = cg->getMass(); |
622 | > | for(atom = cg->beginAtom(ai); atom != NULL; atom = cg->nextAtom(ai)) { |
623 | > | mfact.push_back(atom->getMass()/totalMass); |
624 | > | } |
625 | > | |
626 | > | } |
627 | } | |
628 | < | } |
628 | > | |
629 | > | //fill ident array of local atoms (it is actually ident of AtomType, it is so confusing !!!) |
630 | > | std::vector<int> identArray; |
631 | > | |
632 | > | //to avoid memory reallocation, reserve enough space identArray |
633 | > | identArray.reserve(getNAtoms()); |
634 | ||
635 | < | #ifdef IS_MPI |
636 | < | MPI_Allreduce(&ndfRaw_local,&ndfRaw,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
637 | < | #else |
638 | < | ndfRaw = ndfRaw_local; |
639 | < | #endif |
635 | > | for(mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
636 | > | for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
637 | > | identArray.push_back(atom->getIdent()); |
638 | > | } |
639 | > | } |
640 | ||
641 | < | return ndfRaw; |
642 | < | } |
641 | > | //fill molMembershipArray |
642 | > | //molMembershipArray is filled by SimCreator |
643 | > | std::vector<int> molMembershipArray(nGlobalAtoms_); |
644 | > | for (int i = 0; i < nGlobalAtoms_; i++) { |
645 | > | molMembershipArray[i] = globalMolMembership_[i] + 1; |
646 | > | } |
647 | > | |
648 | > | //setup fortran simulation |
649 | > | //gloalExcludes and molMembershipArray should go away (They are never used) |
650 | > | //why the hell fortran need to know molecule? |
651 | > | //OOPSE = Object-Obfuscated Parallel Simulation Engine |
652 | > | int nGlobalExcludes = 0; |
653 | > | int* globalExcludes = NULL; |
654 | > | int* excludeList = exclude_.getExcludeList(); |
655 | > | setFortranSim( &fInfo_, &nGlobalAtoms_, &nAtoms_, &identArray[0], &nExclude, excludeList , |
656 | > | &nGlobalExcludes, globalExcludes, &molMembershipArray[0], |
657 | > | &mfact[0], &nCutoffGroups_, &fortranGlobalGroupMembership[0], &isError); |
658 | ||
659 | < | int SimInfo::getNDFtranslational() { |
388 | < | int ndfTrans_local; |
659 | > | if( isError ){ |
660 | ||
661 | < | ndfTrans_local = 3 * integrableObjects.size() - n_constraints; |
661 | > | sprintf( painCave.errMsg, |
662 | > | "There was an error setting the simulation information in fortran.\n" ); |
663 | > | painCave.isFatal = 1; |
664 | > | painCave.severity = OOPSE_ERROR; |
665 | > | simError(); |
666 | > | } |
667 | ||
668 | + | #ifdef IS_MPI |
669 | + | sprintf( checkPointMsg, |
670 | + | "succesfully sent the simulation information to fortran.\n"); |
671 | + | MPIcheckPoint(); |
672 | + | #endif // is_mpi |
673 | + | } |
674 | ||
675 | + | |
676 | #ifdef IS_MPI | |
677 | < | MPI_Allreduce(&ndfTrans_local,&ndfTrans,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
678 | < | #else |
679 | < | ndfTrans = ndfTrans_local; |
680 | < | #endif |
677 | > | void SimInfo::setupFortranParallel() { |
678 | > | |
679 | > | //SimInfo is responsible for creating localToGlobalAtomIndex and localToGlobalGroupIndex |
680 | > | std::vector<int> localToGlobalAtomIndex(getNAtoms(), 0); |
681 | > | std::vector<int> localToGlobalCutoffGroupIndex; |
682 | > | SimInfo::MoleculeIterator mi; |
683 | > | Molecule::AtomIterator ai; |
684 | > | Molecule::CutoffGroupIterator ci; |
685 | > | Molecule* mol; |
686 | > | Atom* atom; |
687 | > | CutoffGroup* cg; |
688 | > | mpiSimData parallelData; |
689 | > | int isError; |
690 | ||
691 | < | ndfTrans = ndfTrans - 3 - nZconstraints; |
691 | > | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
692 | ||
693 | < | return ndfTrans; |
694 | < | } |
693 | > | //local index(index in DataStorge) of atom is important |
694 | > | for (atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) { |
695 | > | localToGlobalAtomIndex[atom->getLocalIndex()] = atom->getGlobalIndex() + 1; |
696 | > | } |
697 | ||
698 | < | int SimInfo::getTotIntegrableObjects() { |
699 | < | int nObjs_local; |
700 | < | int nObjs; |
698 | > | //local index of cutoff group is trivial, it only depends on the order of travesing |
699 | > | for (cg = mol->beginCutoffGroup(ci); cg != NULL; cg = mol->nextCutoffGroup(ci)) { |
700 | > | localToGlobalCutoffGroupIndex.push_back(cg->getGlobalIndex() + 1); |
701 | > | } |
702 | > | |
703 | > | } |
704 | ||
705 | < | nObjs_local = integrableObjects.size(); |
705 | > | //fill up mpiSimData struct |
706 | > | parallelData.nMolGlobal = getNGlobalMolecules(); |
707 | > | parallelData.nMolLocal = getNMolecules(); |
708 | > | parallelData.nAtomsGlobal = getNGlobalAtoms(); |
709 | > | parallelData.nAtomsLocal = getNAtoms(); |
710 | > | parallelData.nGroupsGlobal = getNGlobalCutoffGroups(); |
711 | > | parallelData.nGroupsLocal = getNCutoffGroups(); |
712 | > | parallelData.myNode = worldRank; |
713 | > | MPI_Comm_size(MPI_COMM_WORLD, &(parallelData.nProcessors)); |
714 | ||
715 | + | //pass mpiSimData struct and index arrays to fortran |
716 | + | setFsimParallel(¶llelData, &(parallelData.nAtomsLocal), |
717 | + | &localToGlobalAtomIndex[0], &(parallelData.nGroupsLocal), |
718 | + | &localToGlobalCutoffGroupIndex[0], &isError); |
719 | ||
720 | < | #ifdef IS_MPI |
721 | < | MPI_Allreduce(&nObjs_local,&nObjs,1,MPI_INT,MPI_SUM, MPI_COMM_WORLD); |
722 | < | #else |
723 | < | nObjs = nObjs_local; |
724 | < | #endif |
720 | > | if (isError) { |
721 | > | sprintf(painCave.errMsg, |
722 | > | "mpiRefresh errror: fortran didn't like something we gave it.\n"); |
723 | > | painCave.isFatal = 1; |
724 | > | simError(); |
725 | > | } |
726 | ||
727 | + | sprintf(checkPointMsg, " mpiRefresh successful.\n"); |
728 | + | MPIcheckPoint(); |
729 | ||
730 | < | return nObjs; |
730 | > | |
731 | } | |
732 | ||
733 | < | void SimInfo::refreshSim(){ |
733 | > | #endif |
734 | ||
735 | < | simtype fInfo; |
424 | < | int isError; |
425 | < | int n_global; |
426 | < | int* excl; |
735 | > | double SimInfo::calcMaxCutoffRadius() { |
736 | ||
428 | – | fInfo.dielect = 0.0; |
737 | ||
738 | < | if( useDipoles ){ |
739 | < | if( useReactionField )fInfo.dielect = dielectric; |
740 | < | } |
738 | > | std::set<AtomType*> atomTypes; |
739 | > | std::set<AtomType*>::iterator i; |
740 | > | std::vector<double> cutoffRadius; |
741 | ||
742 | < | fInfo.SIM_uses_PBC = usePBC; |
743 | < | //fInfo.SIM_uses_LJ = 0; |
436 | < | fInfo.SIM_uses_LJ = useLJ; |
437 | < | fInfo.SIM_uses_sticky = useSticky; |
438 | < | //fInfo.SIM_uses_sticky = 0; |
439 | < | fInfo.SIM_uses_charges = useCharges; |
440 | < | fInfo.SIM_uses_dipoles = useDipoles; |
441 | < | //fInfo.SIM_uses_dipoles = 0; |
442 | < | fInfo.SIM_uses_RF = useReactionField; |
443 | < | //fInfo.SIM_uses_RF = 0; |
444 | < | fInfo.SIM_uses_GB = useGB; |
445 | < | fInfo.SIM_uses_EAM = useEAM; |
742 | > | //get the unique atom types |
743 | > | atomTypes = getUniqueAtomTypes(); |
744 | ||
745 | < | n_exclude = excludes->getSize(); |
746 | < | excl = excludes->getFortranArray(); |
747 | < | |
745 | > | //query the max cutoff radius among these atom types |
746 | > | for (i = atomTypes.begin(); i != atomTypes.end(); ++i) { |
747 | > | cutoffRadius.push_back(forceField_->getRcutFromAtomType(*i)); |
748 | > | } |
749 | > | |
750 | > | double maxCutoffRadius = *(std::max_element(cutoffRadius.begin(), cutoffRadius.end())); |
751 | #ifdef IS_MPI | |
752 | < | n_global = mpiSim->getNAtomsGlobal(); |
452 | < | #else |
453 | < | n_global = n_atoms; |
752 | > | //pick the max cutoff radius among the processors |
753 | #endif | |
455 | – | |
456 | – | isError = 0; |
457 | – | |
458 | – | getFortranGroupArrays(this, FglobalGroupMembership, mfact); |
459 | – | //it may not be a good idea to pass the address of first element in vector |
460 | – | //since c++ standard does not require vector to be stored continuously in meomory |
461 | – | //Most of the compilers will organize the memory of vector continuously |
462 | – | setFsimulation( &fInfo, &n_global, &n_atoms, identArray, &n_exclude, excl, |
463 | – | &nGlobalExcludes, globalExcludes, molMembershipArray, |
464 | – | &mfact[0], &ngroup, &FglobalGroupMembership[0], &isError); |
754 | ||
755 | < | if( isError ){ |
467 | < | |
468 | < | sprintf( painCave.errMsg, |
469 | < | "There was an error setting the simulation information in fortran.\n" ); |
470 | < | painCave.isFatal = 1; |
471 | < | painCave.severity = OOPSE_ERROR; |
472 | < | simError(); |
473 | < | } |
474 | < | |
475 | < | #ifdef IS_MPI |
476 | < | sprintf( checkPointMsg, |
477 | < | "succesfully sent the simulation information to fortran.\n"); |
478 | < | MPIcheckPoint(); |
479 | < | #endif // is_mpi |
480 | < | |
481 | < | this->ndf = this->getNDF(); |
482 | < | this->ndfRaw = this->getNDFraw(); |
483 | < | this->ndfTrans = this->getNDFtranslational(); |
755 | > | return maxCutoffRadius; |
756 | } | |
757 | ||
758 | < | void SimInfo::setDefaultRcut( double theRcut ){ |
759 | < | |
760 | < | haveRcut = 1; |
761 | < | rCut = theRcut; |
762 | < | rList = rCut + 1.0; |
763 | < | |
764 | < | notifyFortranCutOffs( &rCut, &rSw, &rList ); |
758 | > | void SimInfo::getCutoff(double& rcut, double& rsw) { |
759 | > | |
760 | > | if (fInfo_.SIM_uses_Charges | fInfo_.SIM_uses_Dipoles | fInfo_.SIM_uses_RF) { |
761 | > | |
762 | > | if (!simParams_->haveRcut()){ |
763 | > | sprintf(painCave.errMsg, |
764 | > | "SimCreator Warning: No value was set for the cutoffRadius.\n" |
765 | > | "\tOOPSE will use a default value of 15.0 angstroms" |
766 | > | "\tfor the cutoffRadius.\n"); |
767 | > | painCave.isFatal = 0; |
768 | > | simError(); |
769 | > | rcut = 15.0; |
770 | > | } else{ |
771 | > | rcut = simParams_->getRcut(); |
772 | > | } |
773 | > | |
774 | > | if (!simParams_->haveRsw()){ |
775 | > | sprintf(painCave.errMsg, |
776 | > | "SimCreator Warning: No value was set for switchingRadius.\n" |
777 | > | "\tOOPSE will use a default value of\n" |
778 | > | "\t0.95 * cutoffRadius for the switchingRadius\n"); |
779 | > | painCave.isFatal = 0; |
780 | > | simError(); |
781 | > | rsw = 0.95 * rcut; |
782 | > | } else{ |
783 | > | rsw = simParams_->getRsw(); |
784 | > | } |
785 | > | |
786 | > | } else { |
787 | > | // if charge, dipole or reaction field is not used and the cutofff radius is not specified in |
788 | > | //meta-data file, the maximum cutoff radius calculated from forcefiled will be used |
789 | > | |
790 | > | if (simParams_->haveRcut()) { |
791 | > | rcut = simParams_->getRcut(); |
792 | > | } else { |
793 | > | //set cutoff radius to the maximum cutoff radius based on atom types in the whole system |
794 | > | rcut = calcMaxCutoffRadius(); |
795 | > | } |
796 | > | |
797 | > | if (simParams_->haveRsw()) { |
798 | > | rsw = simParams_->getRsw(); |
799 | > | } else { |
800 | > | rsw = rcut; |
801 | > | } |
802 | > | |
803 | > | } |
804 | } | |
805 | ||
806 | < | void SimInfo::setDefaultRcut( double theRcut, double theRsw ){ |
806 | > | void SimInfo::setupCutoff() { |
807 | > | getCutoff(rcut_, rsw_); |
808 | > | double rnblist = rcut_ + 1; // skin of neighbor list |
809 | ||
810 | < | rSw = theRsw; |
811 | < | setDefaultRcut( theRcut ); |
810 | > | //Pass these cutoff radius etc. to fortran. This function should be called once and only once |
811 | > | notifyFortranCutoffs(&rcut_, &rsw_, &rnblist); |
812 | } | |
813 | ||
814 | + | void SimInfo::addProperty(GenericData* genData) { |
815 | + | properties_.addProperty(genData); |
816 | + | } |
817 | ||
818 | < | void SimInfo::checkCutOffs( void ){ |
819 | < | |
504 | < | if( boxIsInit ){ |
505 | < | |
506 | < | //we need to check cutOffs against the box |
507 | < | |
508 | < | if( rCut > maxCutoff ){ |
509 | < | sprintf( painCave.errMsg, |
510 | < | "cutoffRadius is too large for the current periodic box.\n" |
511 | < | "\tCurrent Value of cutoffRadius = %G at time %G\n " |
512 | < | "\tThis is larger than half of at least one of the\n" |
513 | < | "\tperiodic box vectors. Right now, the Box matrix is:\n" |
514 | < | "\n" |
515 | < | "\t[ %G %G %G ]\n" |
516 | < | "\t[ %G %G %G ]\n" |
517 | < | "\t[ %G %G %G ]\n", |
518 | < | rCut, currentTime, |
519 | < | Hmat[0][0], Hmat[0][1], Hmat[0][2], |
520 | < | Hmat[1][0], Hmat[1][1], Hmat[1][2], |
521 | < | Hmat[2][0], Hmat[2][1], Hmat[2][2]); |
522 | < | painCave.severity = OOPSE_ERROR; |
523 | < | painCave.isFatal = 1; |
524 | < | simError(); |
525 | < | } |
526 | < | } else { |
527 | < | // initialize this stuff before using it, OK? |
528 | < | sprintf( painCave.errMsg, |
529 | < | "Trying to check cutoffs without a box.\n" |
530 | < | "\tOOPSE should have better programmers than that.\n" ); |
531 | < | painCave.severity = OOPSE_ERROR; |
532 | < | painCave.isFatal = 1; |
533 | < | simError(); |
534 | < | } |
535 | < | |
818 | > | void SimInfo::removeProperty(const std::string& propName) { |
819 | > | properties_.removeProperty(propName); |
820 | } | |
821 | ||
822 | < | void SimInfo::addProperty(GenericData* prop){ |
822 | > | void SimInfo::clearProperties() { |
823 | > | properties_.clearProperties(); |
824 | > | } |
825 | ||
826 | < | map<string, GenericData*>::iterator result; |
827 | < | result = properties.find(prop->getID()); |
828 | < | |
543 | < | //we can't simply use properties[prop->getID()] = prop, |
544 | < | //it will cause memory leak if we already contain a propery which has the same name of prop |
545 | < | |
546 | < | if(result != properties.end()){ |
547 | < | |
548 | < | delete (*result).second; |
549 | < | (*result).second = prop; |
826 | > | std::vector<std::string> SimInfo::getPropertyNames() { |
827 | > | return properties_.getPropertyNames(); |
828 | > | } |
829 | ||
830 | < | } |
831 | < | else{ |
830 | > | std::vector<GenericData*> SimInfo::getProperties() { |
831 | > | return properties_.getProperties(); |
832 | > | } |
833 | ||
834 | < | properties[prop->getID()] = prop; |
834 | > | GenericData* SimInfo::getPropertyByName(const std::string& propName) { |
835 | > | return properties_.getPropertyByName(propName); |
836 | > | } |
837 | ||
838 | < | } |
838 | > | void SimInfo::setSnapshotManager(SnapshotManager* sman) { |
839 | > | //if (sman_ == sman_) { |
840 | > | // return; |
841 | > | //} |
842 | ||
843 | + | //delete sman_; |
844 | + | sman_ = sman; |
845 | + | |
846 | + | Molecule* mol; |
847 | + | RigidBody* rb; |
848 | + | Atom* atom; |
849 | + | SimInfo::MoleculeIterator mi; |
850 | + | Molecule::RigidBodyIterator rbIter; |
851 | + | Molecule::AtomIterator atomIter;; |
852 | + | |
853 | + | for (mol = beginMolecule(mi); mol != NULL; mol = nextMolecule(mi)) { |
854 | + | |
855 | + | for (atom = mol->beginAtom(atomIter); atom != NULL; atom = mol->nextAtom(atomIter)) { |
856 | + | atom->setSnapshotManager(sman_); |
857 | + | } |
858 | + | |
859 | + | for (rb = mol->beginRigidBody(rbIter); rb != NULL; rb = mol->nextRigidBody(rbIter)) { |
860 | + | rb->setSnapshotManager(sman_); |
861 | + | } |
862 | + | } |
863 | + | |
864 | } | |
865 | ||
866 | < | GenericData* SimInfo::getProperty(const string& propName){ |
866 | > | Vector3d SimInfo::getComVel(){ |
867 | > | SimInfo::MoleculeIterator i; |
868 | > | Molecule* mol; |
869 | > | |
870 | > | Vector3d comVel(0.0); |
871 | > | double totalMass = 0.0; |
872 | > | |
873 | ||
874 | < | map<string, GenericData*>::iterator result; |
875 | < | |
876 | < | //string lowerCaseName = (); |
877 | < | |
878 | < | result = properties.find(propName); |
879 | < | |
880 | < | if(result != properties.end()) |
881 | < | return (*result).second; |
882 | < | else |
883 | < | return NULL; |
874 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
875 | > | double mass = mol->getMass(); |
876 | > | totalMass += mass; |
877 | > | comVel += mass * mol->getComVel(); |
878 | > | } |
879 | > | |
880 | > | #ifdef IS_MPI |
881 | > | double tmpMass = totalMass; |
882 | > | Vector3d tmpComVel(comVel); |
883 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
884 | > | MPI_Allreduce(tmpComVel.getArrayPointer(), comVel.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
885 | > | #endif |
886 | > | |
887 | > | comVel /= totalMass; |
888 | > | |
889 | > | return comVel; |
890 | } | |
891 | ||
892 | + | Vector3d SimInfo::getCom(){ |
893 | + | SimInfo::MoleculeIterator i; |
894 | + | Molecule* mol; |
895 | ||
896 | < | void SimInfo::getFortranGroupArrays(SimInfo* info, |
897 | < | vector<int>& FglobalGroupMembership, |
898 | < | vector<double>& mfact){ |
899 | < | |
900 | < | Molecule* myMols; |
901 | < | Atom** myAtoms; |
902 | < | int numAtom; |
903 | < | double mtot; |
583 | < | int numMol; |
584 | < | int numCutoffGroups; |
585 | < | CutoffGroup* myCutoffGroup; |
586 | < | vector<CutoffGroup*>::iterator iterCutoff; |
587 | < | Atom* cutoffAtom; |
588 | < | vector<Atom*>::iterator iterAtom; |
589 | < | int atomIndex; |
590 | < | double totalMass; |
591 | < | |
592 | < | mfact.clear(); |
593 | < | FglobalGroupMembership.clear(); |
594 | < | |
896 | > | Vector3d com(0.0); |
897 | > | double totalMass = 0.0; |
898 | > | |
899 | > | for (mol = beginMolecule(i); mol != NULL; mol = nextMolecule(i)) { |
900 | > | double mass = mol->getMass(); |
901 | > | totalMass += mass; |
902 | > | com += mass * mol->getCom(); |
903 | > | } |
904 | ||
596 | – | // Fix the silly fortran indexing problem |
905 | #ifdef IS_MPI | |
906 | < | numAtom = mpiSim->getNAtomsGlobal(); |
907 | < | #else |
908 | < | numAtom = n_atoms; |
906 | > | double tmpMass = totalMass; |
907 | > | Vector3d tmpCom(com); |
908 | > | MPI_Allreduce(&tmpMass,&totalMass,1,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
909 | > | MPI_Allreduce(tmpCom.getArrayPointer(), com.getArrayPointer(),3,MPI_DOUBLE,MPI_SUM, MPI_COMM_WORLD); |
910 | #endif | |
602 | – | for (int i = 0; i < numAtom; i++) |
603 | – | FglobalGroupMembership.push_back(globalGroupMembership[i] + 1); |
604 | – | |
911 | ||
912 | < | myMols = info->molecules; |
607 | < | numMol = info->n_mol; |
608 | < | for(int i = 0; i < numMol; i++){ |
609 | < | numCutoffGroups = myMols[i].getNCutoffGroups(); |
610 | < | for(myCutoffGroup =myMols[i].beginCutoffGroup(iterCutoff); |
611 | < | myCutoffGroup != NULL; |
612 | < | myCutoffGroup =myMols[i].nextCutoffGroup(iterCutoff)){ |
912 | > | com /= totalMass; |
913 | ||
914 | < | totalMass = myCutoffGroup->getMass(); |
615 | < | |
616 | < | for(cutoffAtom = myCutoffGroup->beginAtom(iterAtom); |
617 | < | cutoffAtom != NULL; |
618 | < | cutoffAtom = myCutoffGroup->nextAtom(iterAtom)){ |
619 | < | mfact.push_back(cutoffAtom->getMass()/totalMass); |
620 | < | } |
621 | < | } |
622 | < | } |
914 | > | return com; |
915 | ||
916 | + | } |
917 | + | |
918 | + | std::ostream& operator <<(std::ostream& o, SimInfo& info) { |
919 | + | |
920 | + | return o; |
921 | } | |
922 | + | |
923 | + | }//end namespace oopse |
924 | + |
– | Removed lines |
+ | Added lines |
< | Changed lines |
> | Changed lines |