ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/integrators/LDForceManager.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/integrators/LDForceManager.cpp (file contents):
Revision 2634 by tim, Fri Mar 17 23:20:35 2006 UTC vs.
Revision 2733 by gezelter, Tue Apr 25 02:09:01 2006 UTC

# Line 46 | Line 46 | namespace oopse {
46  
47    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
48      Globals* simParams = info->getSimParams();
49 +    
50      std::map<std::string, HydroProp> hydroPropMap;
51      if (simParams->haveHydroPropFile()) {
52 <        hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
52 >      hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
53      } else {
54 <        //error
54 >      sprintf( painCave.errMsg,
55 >               "HydroPropFile must be set if Langevin Dynamics is specified.\n");
56 >      painCave.severity = OOPSE_ERROR;
57 >      painCave.isFatal = 1;
58 >      simError();  
59      }
60 +    
61 +    sphericalBoundaryConditions_ = false;
62 +    if (simParams->getUseSphericalBoundaryConditions()) {
63 +      sphericalBoundaryConditions_ = true;
64 +      if (simParams->haveLangevinBufferRadius()) {
65 +        langevinBufferRadius_ = simParams->getLangevinBufferRadius();
66 +      } else {
67 +        sprintf( painCave.errMsg,
68 +                 "langevinBufferRadius must be specified "
69 +                 "when useSphericalBoundaryConditions is turned on.\n");
70 +        painCave.severity = OOPSE_ERROR;
71 +        painCave.isFatal = 1;
72 +        simError();  
73 +      }
74 +    
75 +      if (simParams->haveFrozenBufferRadius()) {
76 +        frozenBufferRadius_ = simParams->getFrozenBufferRadius();
77 +      } else {
78 +        sprintf( painCave.errMsg,
79 +                 "frozenBufferRadius must be specified "
80 +                 "when useSphericalBoundaryConditions is turned on.\n");
81 +        painCave.severity = OOPSE_ERROR;
82 +        painCave.isFatal = 1;
83 +        simError();  
84 +      }
85  
86 +      if (frozenBufferRadius_ < langevinBufferRadius_) {
87 +        sprintf( painCave.errMsg,
88 +                 "frozenBufferRadius has been set smaller than the "
89 +                 "langevinBufferRadius.  This is probably an error.\n");
90 +        painCave.severity = OOPSE_WARNING;
91 +        painCave.isFatal = 0;
92 +        simError();  
93 +      }
94 +    }
95 +      
96      SimInfo::MoleculeIterator i;
97      Molecule::IntegrableObjectIterator  j;
98      Molecule* mol;
99      StuntDouble* integrableObject;
100      for (mol = info->beginMolecule(i); mol != NULL; mol = info->nextMolecule(i)) {
101        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
102 <              integrableObject = mol->nextIntegrableObject(j)) {
103 <            std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
104 <            if (iter != hydroPropMap.end()) {
105 <                hydroProps_.push_back(iter->second);
106 <            } else {
107 <                //error
108 <            }
109 <            
110 <           }
102 >           integrableObject = mol->nextIntegrableObject(j)) {
103 >        std::map<std::string, HydroProp>::iterator iter = hydroPropMap.find(integrableObject->getType());
104 >        if (iter != hydroPropMap.end()) {
105 >          hydroProps_.push_back(iter->second);
106 >        } else {
107 >          sprintf( painCave.errMsg,
108 >                   "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
109 >          painCave.severity = OOPSE_ERROR;
110 >          painCave.isFatal = 1;
111 >          simError();  
112 >        }
113 >        
114 >      }
115      }
116      variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
117    }
# Line 75 | Line 119 | namespace oopse {
119      std::map<std::string, HydroProp> props;
120      std::ifstream ifs(filename.c_str());
121      if (ifs.is_open()) {
122 <
122 >      
123      }
124 <
124 >    
125      const unsigned int BufferSize = 65535;
126      char buffer[BufferSize];  
83 –    Mat3x3d Ddtt;
84 –    Mat3x3d Ddtr;
85 –    Mat3x3d Ddrr;
127      while (ifs.getline(buffer, BufferSize)) {
128 <        StringTokenizer tokenizer(buffer);
129 <        HydroProp currProp;
130 <        if (tokenizer.countTokens() >= 40) {
131 <            std::string atomName = tokenizer.nextToken();
132 <            currProp.cor[0] = tokenizer.nextTokenAsDouble();
133 <            currProp.cor[1] = tokenizer.nextTokenAsDouble();
134 <            currProp.cor[2] = tokenizer.nextTokenAsDouble();
135 <
136 <
137 <            Ddtt(0,0) = tokenizer.nextTokenAsDouble();
138 <            Ddtt(0,1) = tokenizer.nextTokenAsDouble();
139 <            Ddtt(0,2) = tokenizer.nextTokenAsDouble();
140 <            Ddtt(1,0) = tokenizer.nextTokenAsDouble();
141 <            Ddtt(1,1) = tokenizer.nextTokenAsDouble();
142 <            Ddtt(1,2) = tokenizer.nextTokenAsDouble();
143 <            Ddtt(2,0) = tokenizer.nextTokenAsDouble();
144 <            Ddtt(2,1) = tokenizer.nextTokenAsDouble();
145 <            Ddtt(2,2) = tokenizer.nextTokenAsDouble();
146 <
147 <            Ddtr(0,0) = tokenizer.nextTokenAsDouble();
148 <            Ddtr(0,1) = tokenizer.nextTokenAsDouble();
149 <            Ddtr(0,2) = tokenizer.nextTokenAsDouble();
150 <            Ddtr(1,0) = tokenizer.nextTokenAsDouble();
151 <            Ddtr(1,1) = tokenizer.nextTokenAsDouble();
152 <            Ddtr(1,2) = tokenizer.nextTokenAsDouble();
153 <            Ddtr(2,0) = tokenizer.nextTokenAsDouble();
154 <            Ddtr(2,1) = tokenizer.nextTokenAsDouble();
155 <            Ddtr(2,2) = tokenizer.nextTokenAsDouble();
156 <
157 <            Ddrr(0,0) = tokenizer.nextTokenAsDouble();
158 <            Ddrr(0,1) = tokenizer.nextTokenAsDouble();
159 <            Ddrr(0,2) = tokenizer.nextTokenAsDouble();
160 <            Ddrr(1,0) = tokenizer.nextTokenAsDouble();
161 <            Ddrr(1,1) = tokenizer.nextTokenAsDouble();
162 <            Ddrr(1,2) = tokenizer.nextTokenAsDouble();
163 <            Ddrr(2,0) = tokenizer.nextTokenAsDouble();
164 <            Ddrr(2,1) = tokenizer.nextTokenAsDouble();
165 <            Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
166 <
167 <            currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
168 <            currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
169 <            currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
170 <            currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
171 <            currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
172 <            currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
173 <            currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
174 <            currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
175 <            currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
176 <
177 <            currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
178 <            currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
179 <            currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
180 <            currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
181 <            currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
182 <            currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
183 <            currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
184 <            currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
145 <        
146 <            currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
147 <            currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
148 <            currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
149 <            currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
150 <            currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
151 <            currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
152 <            currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
153 <            currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
154 <            currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
155 <
156 <            currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
157 <            currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
158 <            currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
159 <            currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
160 <            currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
161 <            currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
162 <            currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
163 <            currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
164 <            currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
165 <
166 <            SquareMatrix<double, 6> Xir;
167 <            Xir.setSubMatrix(0, 0, currProp.Xirtt);
168 <            Xir.setSubMatrix(0, 3, currProp.Xirrt);
169 <            Xir.setSubMatrix(3, 0, currProp.Xirtr);
170 <            Xir.setSubMatrix(3, 3, currProp.Xirrr);
171 <            CholeskyDecomposition(Xir, currProp.S);            
172 <
173 <            props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
174 <        }
128 >      StringTokenizer tokenizer(buffer);
129 >      HydroProp currProp;
130 >      if (tokenizer.countTokens() >= 40) {
131 >        std::string atomName = tokenizer.nextToken();
132 >        currProp.cor[0] = tokenizer.nextTokenAsDouble();
133 >        currProp.cor[1] = tokenizer.nextTokenAsDouble();
134 >        currProp.cor[2] = tokenizer.nextTokenAsDouble();
135 >        
136 >        currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
137 >        currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
138 >        currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
139 >        currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
140 >        currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
141 >        currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
142 >        currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
143 >        currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
144 >        currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
145 >        
146 >        currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
147 >        currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
148 >        currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
149 >        currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
150 >        currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
151 >        currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
152 >        currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
153 >        currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
154 >        currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
155 >        
156 >        currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
157 >        currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
158 >        currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
159 >        currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
160 >        currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
161 >        currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
162 >        currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
163 >        currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
164 >        currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
165 >        
166 >        currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
167 >        currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
168 >        currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
169 >        currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
170 >        currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
171 >        currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
172 >        currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
173 >        currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
174 >        currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
175 >        
176 >        SquareMatrix<double, 6> Xir;
177 >        Xir.setSubMatrix(0, 0, currProp.Xirtt);
178 >        Xir.setSubMatrix(0, 3, currProp.Xirrt);
179 >        Xir.setSubMatrix(3, 0, currProp.Xirtr);
180 >        Xir.setSubMatrix(3, 3, currProp.Xirrr);
181 >        CholeskyDecomposition(Xir, currProp.S);            
182 >        
183 >        props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
184 >      }
185      }
186 <
186 >    
187      return props;
188    }
189    
# Line 191 | Line 201 | namespace oopse {
201      Vector3d ji;
202      double mass;
203      unsigned int index = 0;
204 +    bool doLangevinForces;
205 +    bool freezeMolecule;
206 +    int fdf;
207 +    
208 +    fdf = 0;
209      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
210 <      for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
211 <           integrableObject = mol->nextIntegrableObject(j)) {
212 <
210 >      
211 >      if (sphericalBoundaryConditions_) {
212 >        
213 >        Vector3d molPos = mol->getCom();
214 >        double molRad = molPos.length();
215 >        
216 >        doLangevinForces = false;
217 >        freezeMolecule = false;
218 >        
219 >        if (molRad > langevinBufferRadius_) {
220 >          doLangevinForces = true;
221 >          freezeMolecule = false;
222 >        }
223 >        if (molRad > frozenBufferRadius_) {
224 >          doLangevinForces = false;
225 >          freezeMolecule = true;
226 >        }
227 >      }
228 >      
229 >      if (doLangevinForces) {
230 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
231 >             integrableObject = mol->nextIntegrableObject(j)) {
232 >          
233            vel =integrableObject->getVel();
234            if (integrableObject->isDirectional()){
235 <             //calculate angular velocity in lab frame
236 <             Mat3x3d I = integrableObject->getI();
237 <             Vector3d angMom = integrableObject->getJ();
238 <             Vector3d omega;
239 <
240 <             if (integrableObject->isLinear()) {
241 <                int linearAxis = integrableObject->linearAxis();
242 <                int l = (linearAxis +1 )%3;
243 <                int m = (linearAxis +2 )%3;
244 <                omega[l] = angMom[l] /I(l, l);
245 <                omega[m] = angMom[m] /I(m, m);
246 <                
247 <             } else {
248 <                 omega[0] = angMom[0] /I(0, 0);
249 <                 omega[1] = angMom[1] /I(1, 1);
250 <                 omega[2] = angMom[2] /I(2, 2);
251 <             }
252 <
253 <             //apply friction force and torque at center of resistance
254 <             A = integrableObject->getA();
255 <             Atrans = A.transpose();
256 <             Vector3d rcr = Atrans * hydroProps_[index].cor;  
257 <             Vector3d vcdLab = vel + cross(omega, rcr);
258 <             Vector3d vcdBody = A* vcdLab;
259 <             Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
260 <             Vector3d frictionForceLab = Atrans*frictionForceBody;
261 <             integrableObject->addFrc(frictionForceLab);
262 <             Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
263 <             Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
264 <             integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
265 <
266 <             //apply random force and torque at center of resistance
267 <             Vector3d randomForceBody;
268 <             Vector3d randomTorqueBody;
269 <             genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
270 <             Vector3d randomForceLab = Atrans*randomForceBody;
271 <             Vector3d randomTorqueLab = Atrans* randomTorqueBody;
272 <             integrableObject->addFrc(randomForceLab);            
273 <             integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
274 <
235 >            //calculate angular velocity in lab frame
236 >            Mat3x3d I = integrableObject->getI();
237 >            Vector3d angMom = integrableObject->getJ();
238 >            Vector3d omega;
239 >            
240 >            if (integrableObject->isLinear()) {
241 >              int linearAxis = integrableObject->linearAxis();
242 >              int l = (linearAxis +1 )%3;
243 >              int m = (linearAxis +2 )%3;
244 >              omega[l] = angMom[l] /I(l, l);
245 >              omega[m] = angMom[m] /I(m, m);
246 >              
247 >            } else {
248 >              omega[0] = angMom[0] /I(0, 0);
249 >              omega[1] = angMom[1] /I(1, 1);
250 >              omega[2] = angMom[2] /I(2, 2);
251 >            }
252 >            
253 >            //apply friction force and torque at center of resistance
254 >            A = integrableObject->getA();
255 >            Atrans = A.transpose();
256 >            Vector3d rcr = Atrans * hydroProps_[index].cor;  
257 >            Vector3d vcdLab = vel + cross(omega, rcr);
258 >            Vector3d vcdBody = A* vcdLab;
259 >            Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
260 >            Vector3d frictionForceLab = Atrans*frictionForceBody;
261 >            integrableObject->addFrc(frictionForceLab);
262 >            Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
263 >            Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
264 >            integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
265 >            
266 >            //apply random force and torque at center of resistance
267 >            Vector3d randomForceBody;
268 >            Vector3d randomTorqueBody;
269 >            genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
270 >            Vector3d randomForceLab = Atrans*randomForceBody;
271 >            Vector3d randomTorqueLab = Atrans* randomTorqueBody;
272 >            integrableObject->addFrc(randomForceLab);            
273 >            integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
274 >            
275            } else {
276 <             //spheric atom
277 <             Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
278 <             Vector3d randomForce;
279 <             Vector3d randomTorque;
280 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
281 <
282 <             integrableObject->addFrc(frictionForce+randomForce);            
276 >            //spherical atom
277 >            Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
278 >            Vector3d randomForce;
279 >            Vector3d randomTorque;
280 >            genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
281 >            
282 >            integrableObject->addFrc(frictionForce+randomForce);            
283            }
284 <
285 <        ++index;
284 >          
285 >          ++index;
286      
287 +        }
288        }
289 <    }    
290 <
291 <    ForceManager::postCalculation();
292 <
293 <
294 <
289 >      if (freezeMolecule)
290 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
291 >             integrableObject = mol->nextIntegrableObject(j)) {          
292 >          fdf += integrableObject->freeze();
293 >        }
294 >    }
295 >    
296 >    info_->setFdf(fdf);
297 >    
298 >    ForceManager::postCalculation();  
299    }
300  
301   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines