ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/integrators/LDForceManager.cpp
(Generate patch)

Comparing trunk/src/integrators/LDForceManager.cpp (file contents):
Revision 895 by tim, Mon Mar 13 22:42:40 2006 UTC vs.
Revision 909 by tim, Tue Mar 21 00:26:55 2006 UTC

# Line 41 | Line 41
41   #include <fstream>
42   #include "integrators/LDForceManager.hpp"
43   #include "math/CholeskyDecomposition.hpp"
44 + #include "utils/OOPSEConstant.hpp"
45   namespace oopse {
46  
47    LDForceManager::LDForceManager(SimInfo* info) : ForceManager(info){
# Line 49 | Line 50 | namespace oopse {
50      if (simParams->haveHydroPropFile()) {
51          hydroPropMap = parseFrictionFile(simParams->getHydroPropFile());
52      } else {
53 <        //error
53 >                sprintf( painCave.errMsg,
54 >                       "HydroPropFile keyword must be set if Lagevin Dynamics is used\n");
55 >                painCave.severity = OOPSE_ERROR;
56 >                painCave.isFatal = 1;
57 >                simError();  
58      }
59  
60      SimInfo::MoleculeIterator i;
# Line 63 | Line 68 | namespace oopse {
68              if (iter != hydroPropMap.end()) {
69                  hydroProps_.push_back(iter->second);
70              } else {
71 <                //error
71 >                sprintf( painCave.errMsg,
72 >                       "Can not find resistance tensor for atom [%s]\n", integrableObject->getType().c_str());
73 >                painCave.severity = OOPSE_ERROR;
74 >                painCave.isFatal = 1;
75 >                simError();  
76              }
77              
78             }
79      }
80 <    variance_ = 2.0*simParams->getDt();
80 >    variance_ = 2.0 * OOPSEConstant::kb*simParams->getTargetTemp()/simParams->getDt();
81    }
82    std::map<std::string, HydroProp> LDForceManager::parseFrictionFile(const std::string& filename) {
83      std::map<std::string, HydroProp> props;
# Line 82 | Line 91 | namespace oopse {
91      while (ifs.getline(buffer, BufferSize)) {
92          StringTokenizer tokenizer(buffer);
93          HydroProp currProp;
94 <        if (tokenizer.countTokens() >= 67) {
94 >        if (tokenizer.countTokens() >= 40) {
95              std::string atomName = tokenizer.nextToken();
96 <            currProp.cod[0] = tokenizer.nextTokenAsDouble();
97 <            currProp.cod[1] = tokenizer.nextTokenAsDouble();
98 <            currProp.cod[2] = tokenizer.nextTokenAsDouble();
96 >            currProp.cor[0] = tokenizer.nextTokenAsDouble();
97 >            currProp.cor[1] = tokenizer.nextTokenAsDouble();
98 >            currProp.cor[2] = tokenizer.nextTokenAsDouble();
99 >            
100 >            currProp.Xirtt(0,0) = tokenizer.nextTokenAsDouble();
101 >            currProp.Xirtt(0,1) = tokenizer.nextTokenAsDouble();
102 >            currProp.Xirtt(0,2) = tokenizer.nextTokenAsDouble();
103 >            currProp.Xirtt(1,0) = tokenizer.nextTokenAsDouble();
104 >            currProp.Xirtt(1,1) = tokenizer.nextTokenAsDouble();
105 >            currProp.Xirtt(1,2) = tokenizer.nextTokenAsDouble();
106 >            currProp.Xirtt(2,0) = tokenizer.nextTokenAsDouble();
107 >            currProp.Xirtt(2,1) = tokenizer.nextTokenAsDouble();
108 >            currProp.Xirtt(2,2) = tokenizer.nextTokenAsDouble();
109  
110 <            currProp.Ddtt(0,0) = tokenizer.nextTokenAsDouble();
111 <            currProp.Ddtt(0,1) = tokenizer.nextTokenAsDouble();
112 <            currProp.Ddtt(0,2) = tokenizer.nextTokenAsDouble();
113 <            currProp.Ddtt(1,0) = tokenizer.nextTokenAsDouble();
114 <            currProp.Ddtt(1,1) = tokenizer.nextTokenAsDouble();
115 <            currProp.Ddtt(1,2) = tokenizer.nextTokenAsDouble();
116 <            currProp.Ddtt(2,0) = tokenizer.nextTokenAsDouble();
117 <            currProp.Ddtt(2,1) = tokenizer.nextTokenAsDouble();
118 <            currProp.Ddtt(2,2) = tokenizer.nextTokenAsDouble();
110 >            currProp.Xirrt(0,0) = tokenizer.nextTokenAsDouble();
111 >            currProp.Xirrt(0,1) = tokenizer.nextTokenAsDouble();
112 >            currProp.Xirrt(0,2) = tokenizer.nextTokenAsDouble();
113 >            currProp.Xirrt(1,0) = tokenizer.nextTokenAsDouble();
114 >            currProp.Xirrt(1,1) = tokenizer.nextTokenAsDouble();
115 >            currProp.Xirrt(1,2) = tokenizer.nextTokenAsDouble();
116 >            currProp.Xirrt(2,0) = tokenizer.nextTokenAsDouble();
117 >            currProp.Xirrt(2,1) = tokenizer.nextTokenAsDouble();
118 >            currProp.Xirrt(2,2) = tokenizer.nextTokenAsDouble();
119 >        
120 >            currProp.Xirtr(0,0) = tokenizer.nextTokenAsDouble();
121 >            currProp.Xirtr(0,1) = tokenizer.nextTokenAsDouble();
122 >            currProp.Xirtr(0,2) = tokenizer.nextTokenAsDouble();
123 >            currProp.Xirtr(1,0) = tokenizer.nextTokenAsDouble();
124 >            currProp.Xirtr(1,1) = tokenizer.nextTokenAsDouble();
125 >            currProp.Xirtr(1,2) = tokenizer.nextTokenAsDouble();
126 >            currProp.Xirtr(2,0) = tokenizer.nextTokenAsDouble();
127 >            currProp.Xirtr(2,1) = tokenizer.nextTokenAsDouble();
128 >            currProp.Xirtr(2,2) = tokenizer.nextTokenAsDouble();
129  
130 <            currProp.Ddtr(0,0) = tokenizer.nextTokenAsDouble();
131 <            currProp.Ddtr(0,1) = tokenizer.nextTokenAsDouble();
132 <            currProp.Ddtr(0,2) = tokenizer.nextTokenAsDouble();
133 <            currProp.Ddtr(1,0) = tokenizer.nextTokenAsDouble();
134 <            currProp.Ddtr(1,1) = tokenizer.nextTokenAsDouble();
135 <            currProp.Ddtr(1,2) = tokenizer.nextTokenAsDouble();
136 <            currProp.Ddtr(2,0) = tokenizer.nextTokenAsDouble();
137 <            currProp.Ddtr(2,1) = tokenizer.nextTokenAsDouble();
138 <            currProp.Ddtr(2,2) = tokenizer.nextTokenAsDouble();
110 <
111 <            currProp.Ddrr(0,0) = tokenizer.nextTokenAsDouble();
112 <            currProp.Ddrr(0,1) = tokenizer.nextTokenAsDouble();
113 <            currProp.Ddrr(0,2) = tokenizer.nextTokenAsDouble();
114 <            currProp.Ddrr(1,0) = tokenizer.nextTokenAsDouble();
115 <            currProp.Ddrr(1,1) = tokenizer.nextTokenAsDouble();
116 <            currProp.Ddrr(1,2) = tokenizer.nextTokenAsDouble();
117 <            currProp.Ddrr(2,0) = tokenizer.nextTokenAsDouble();
118 <            currProp.Ddrr(2,1) = tokenizer.nextTokenAsDouble();
119 <            currProp.Ddrr(2,2) = tokenizer.nextTokenAsDouble();                
120 <
121 <            currProp.Xidtt(0,0) = tokenizer.nextTokenAsDouble();
122 <            currProp.Xidtt(0,1) = tokenizer.nextTokenAsDouble();
123 <            currProp.Xidtt(0,2) = tokenizer.nextTokenAsDouble();
124 <            currProp.Xidtt(1,0) = tokenizer.nextTokenAsDouble();
125 <            currProp.Xidtt(1,1) = tokenizer.nextTokenAsDouble();
126 <            currProp.Xidtt(1,2) = tokenizer.nextTokenAsDouble();
127 <            currProp.Xidtt(2,0) = tokenizer.nextTokenAsDouble();
128 <            currProp.Xidtt(2,1) = tokenizer.nextTokenAsDouble();
129 <            currProp.Xidtt(2,2) = tokenizer.nextTokenAsDouble();
130 >            currProp.Xirrr(0,0) = tokenizer.nextTokenAsDouble();
131 >            currProp.Xirrr(0,1) = tokenizer.nextTokenAsDouble();
132 >            currProp.Xirrr(0,2) = tokenizer.nextTokenAsDouble();
133 >            currProp.Xirrr(1,0) = tokenizer.nextTokenAsDouble();
134 >            currProp.Xirrr(1,1) = tokenizer.nextTokenAsDouble();
135 >            currProp.Xirrr(1,2) = tokenizer.nextTokenAsDouble();
136 >            currProp.Xirrr(2,0) = tokenizer.nextTokenAsDouble();
137 >            currProp.Xirrr(2,1) = tokenizer.nextTokenAsDouble();
138 >            currProp.Xirrr(2,2) = tokenizer.nextTokenAsDouble();
139  
140 <            currProp.Xidrt(0,0) = tokenizer.nextTokenAsDouble();
141 <            currProp.Xidrt(0,1) = tokenizer.nextTokenAsDouble();
142 <            currProp.Xidrt(0,2) = tokenizer.nextTokenAsDouble();
143 <            currProp.Xidrt(1,0) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xidrt(1,1) = tokenizer.nextTokenAsDouble();
145 <            currProp.Xidrt(1,2) = tokenizer.nextTokenAsDouble();
137 <            currProp.Xidrt(2,0) = tokenizer.nextTokenAsDouble();
138 <            currProp.Xidrt(2,1) = tokenizer.nextTokenAsDouble();
139 <            currProp.Xidrt(2,2) = tokenizer.nextTokenAsDouble();
140 <            
141 <            currProp.Xidtr(0,0) = tokenizer.nextTokenAsDouble();
142 <            currProp.Xidtr(0,1) = tokenizer.nextTokenAsDouble();
143 <            currProp.Xidtr(0,2) = tokenizer.nextTokenAsDouble();
144 <            currProp.Xidtr(1,0) = tokenizer.nextTokenAsDouble();
145 <            currProp.Xidtr(1,1) = tokenizer.nextTokenAsDouble();
146 <            currProp.Xidtr(1,2) = tokenizer.nextTokenAsDouble();
147 <            currProp.Xidtr(2,0) = tokenizer.nextTokenAsDouble();
148 <            currProp.Xidtr(2,1) = tokenizer.nextTokenAsDouble();
149 <            currProp.Xidtr(2,2) = tokenizer.nextTokenAsDouble();
140 >            SquareMatrix<double, 6> Xir;
141 >            Xir.setSubMatrix(0, 0, currProp.Xirtt);
142 >            Xir.setSubMatrix(0, 3, currProp.Xirrt);
143 >            Xir.setSubMatrix(3, 0, currProp.Xirtr);
144 >            Xir.setSubMatrix(3, 3, currProp.Xirrr);
145 >            CholeskyDecomposition(Xir, currProp.S);            
146  
151 –            currProp.Xidrr(0,0) = tokenizer.nextTokenAsDouble();
152 –            currProp.Xidrr(0,1) = tokenizer.nextTokenAsDouble();
153 –            currProp.Xidrr(0,2) = tokenizer.nextTokenAsDouble();
154 –            currProp.Xidrr(1,0) = tokenizer.nextTokenAsDouble();
155 –            currProp.Xidrr(1,1) = tokenizer.nextTokenAsDouble();
156 –            currProp.Xidrr(1,2) = tokenizer.nextTokenAsDouble();
157 –            currProp.Xidrr(2,0) = tokenizer.nextTokenAsDouble();
158 –            currProp.Xidrr(2,1) = tokenizer.nextTokenAsDouble();
159 –            currProp.Xidrr(2,2) = tokenizer.nextTokenAsDouble();
147              props.insert(std::map<std::string, HydroProp>::value_type(atomName, currProp));
148          }
149      }
# Line 173 | Line 160 | namespace oopse {
160      Vector3d pos;
161      Vector3d frc;
162      Mat3x3d A;
163 +    Mat3x3d Atrans;
164      Vector3d Tb;
165      Vector3d ji;
166      double mass;
# Line 201 | Line 189 | namespace oopse {
189                   omega[2] = angMom[2] /I(2, 2);
190               }
191  
192 <             //apply friction force and torque at center of diffusion
192 >             //apply friction force and torque at center of resistance
193               A = integrableObject->getA();
194 <             Vector3d rcd = A.transpose() * hydroProps_[index].cod;  
195 <             Vector3d vcd = vel + cross(omega, rcd);
196 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt * vcd + hydroProps_[index].Xidrt * omega);
197 <             integrableObject->addFrc(frictionForce);
198 <             Vector3d frictionTorque = - (hydroProps_[index].Xidtr * vcd + hydroProps_[index].Xidrr * omega);
199 <             integrableObject->addTrq(frictionTorque);
200 <            
201 <             //apply random force and torque at center of diffustion
202 <             Vector3d randomForce;
203 <             Vector3d randomTorque;
204 <             genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
205 <             integrableObject->addFrc(randomForce);
206 <             integrableObject->addTrq(randomTorque + cross(rcd, randomForce ));
207 <            
194 >             Atrans = A.transpose();
195 >             Vector3d rcr = Atrans * hydroProps_[index].cor;  
196 >             Vector3d vcdLab = vel + cross(omega, rcr);
197 >             Vector3d vcdBody = A* vcdLab;
198 >             Vector3d frictionForceBody = -(hydroProps_[index].Xirtt * vcdBody + hydroProps_[index].Xirrt * omega);
199 >             Vector3d frictionForceLab = Atrans*frictionForceBody;
200 >             integrableObject->addFrc(frictionForceLab);
201 >             Vector3d frictionTorqueBody = - (hydroProps_[index].Xirtr * vcdBody + hydroProps_[index].Xirrr * omega);
202 >             Vector3d frictionTorqueLab = Atrans*frictionTorqueBody;
203 >             integrableObject->addTrq(frictionTorqueLab+ cross(rcr, frictionForceLab));
204 >
205 >             //apply random force and torque at center of resistance
206 >             Vector3d randomForceBody;
207 >             Vector3d randomTorqueBody;
208 >             genRandomForceAndTorque(randomForceBody, randomTorqueBody, index, variance_);
209 >             Vector3d randomForceLab = Atrans*randomForceBody;
210 >             Vector3d randomTorqueLab = Atrans* randomTorqueBody;
211 >             integrableObject->addFrc(randomForceLab);            
212 >             integrableObject->addTrq(randomTorqueLab + cross(rcr, randomForceLab ));            
213 >
214            } else {
215               //spheric atom
216 <             Vector3d frictionForce = -(hydroProps_[index].Xidtt *vel);    
216 >             Vector3d frictionForce = -(hydroProps_[index].Xirtt *vel);    
217               Vector3d randomForce;
218               Vector3d randomTorque;
219               genRandomForceAndTorque(randomForce, randomTorque, index, variance_);
220 +
221               integrableObject->addFrc(frictionForce+randomForce);            
222            }
223  
# Line 238 | Line 233 | void LDForceManager::genRandomForceAndTorque(Vector3d&
233    }
234  
235   void LDForceManager::genRandomForceAndTorque(Vector3d& force, Vector3d& torque, unsigned int index, double variance) {
236 <    SquareMatrix<double, 6> Dd;
237 <    SquareMatrix<double, 6> S;
236 >
237 >
238      Vector<double, 6> Z;
239      Vector<double, 6> generalForce;
240 <    Dd.setSubMatrix(0, 0, hydroProps_[index].Ddtt);
241 <    Dd.setSubMatrix(0, 3, hydroProps_[index].Ddtr.transpose());
247 <    Dd.setSubMatrix(3, 0, hydroProps_[index].Ddtr);
248 <    Dd.setSubMatrix(3, 3, hydroProps_[index].Ddrr);
249 <    CholeskyDecomposition(Dd, S);
250 <    
240 >
241 >        
242      Z[0] = randNumGen_.randNorm(0, variance);
243      Z[1] = randNumGen_.randNorm(0, variance);
244      Z[2] = randNumGen_.randNorm(0, variance);
245      Z[3] = randNumGen_.randNorm(0, variance);
246      Z[4] = randNumGen_.randNorm(0, variance);
247      Z[5] = randNumGen_.randNorm(0, variance);
248 <        
249 <    generalForce = S*Z;
248 >    
249 >
250 >    generalForce = hydroProps_[index].S*Z;
251 >    
252      force[0] = generalForce[0];
253      force[1] = generalForce[1];
254      force[2] = generalForce[2];

Diff Legend

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