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