ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/UseTheForce/DUFF.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/UseTheForce/DUFF.cpp (file contents):
Revision 1757 by tim, Fri Nov 19 00:11:33 2004 UTC vs.
Revision 1760 by tim, Fri Nov 19 20:23:26 2004 UTC

# Line 25 | Line 25
25  
26   #include "UseTheForce/DUFF.hpp"
27   #include "UseTheForce/DarkSide/lj_interface.h"
28 + #include "UseTheForce/DarkSide/charge_interface.h"
29 + #include "UseTheForce/DarkSide/dipole_interface.h"
30 + #include "UseTheForce/DarkSide/sticky_interface.h"
31  
32   namespace oopse {
33  
# Line 33 | Line 36 | ForceField* createDUFF() {
36      return new DUFF();
37   }
38  
39 < //register createDUFF to ForceFieldCreator
39 > //register createDUFF to ForceFieldFactory
40   ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
41  
42  
# Line 44 | Line 47 | ParseState DUFF::getSection(const std::string& section
47          case "AtomTypes" :
48              result = DUFF::AtomTypeSection;
49              break;
50 +        case "DirectionalAtomTypes" :
51 +            result = DUFF::DirectionalAtomTypeSection;
52 +            break;
53  
54          case "BondTypes" :
55              result = DUFF::BondTypeSection;
# Line 70 | Line 76 | void DUFF::parse(const std::string& filename) {
76      std::string line;
77      char buffer[bufferSize];
78      int lineNo = 0;
79 +    int atomIdent = getNAtomType() + 1;  //atom's indent begins from 1 (since only fortran's array begins from 1)
80      ParseState currentSection = DUFF::UnknownSection;
81      
82      while(ffStream.getline(buffer, bufferSize)){
# Line 83 | Line 90 | void DUFF::parse(const std::string& filename) {
90  
91              switch(currentSection) {
92                  case DUFF::AtomTypeSection :
93 <                    parseAtomType(line, lineNo);
93 >                    parseAtomType(line, lineNo, atomIdent);
94                      break;
95  
96 +                case DUFF::DirectionalAtomTypeSection :
97 +                    parseDirectionalAtomType(line, lineNo);
98 +                    break;
99 +                    
100                  case DUFF::BondTypeSection :
101                      parseBondType(line, lineNo);
102                      break;
# Line 136 | Line 147 | void DUFF::parseAtomType(const std::string& line, int
147   }
148  
149  
150 < void DUFF::parseAtomType(const std::string& line, int lineNo){
150 > void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
151      StringTokenizer tokenizer(line);
152 <    
142 <    if (tokenizer.countTokens() >= 4) {
143 <        atomTypeName = tokenizer.nextToken();
144 <        mass = tokenizer.nextTokenAsFloat();
145 <        epsilon = tokenizer.nextTokenAsFloat();
146 <        sigma = tokenizer.nextTokenAsFloat();
152 >    int nTokens = tokenizer.countTokens();    
153  
154 <        atomType = new AtomType();
154 >    //in AtomTypeSection, a line at least contains 5 tokens
155 >    //atomTypeName, is Directional, isLJ, isCharge and mass
156 >    if (nTokens < 5)  {
157 >                      
158 >    } else {
159 >
160 >        std::string atomTypeName = tokenizer.nextToken();
161 >        bool isDirectional = tokenizer.nextTokenAsBool();
162 >        bool isLJ = tokenizer.nextTokenAsBool();
163 >        bool isCharge = tokenizer.nextTokenAsBool();        
164 >        double mass = tokenizer.nextTokenAsDouble();
165 >        double epsilon;
166 >        double sigma;
167 >        double charge;
168 >        nTokens -= 5;
169 >
170 >        //parse epsilon and sigma
171 >        if (isLJ) {
172 >            if (nTokens >= 2) {
173 >                epsilon = tokenizer.nextTokenAsDouble();
174 >                sigma = tokenizer.nextTokenAsDouble();
175 >                nTokens -= 2;
176 >            } else {
177 >
178 >            }
179 >        }
180 >
181 >        //parse charge
182 >        if (isCharge) {
183 >            if (nTokens >= 1) {
184 >                charge = tokenizer.nextTokenAsDouble();
185 >                nTokens -= 1;
186 >            } else {
187 >
188 >            }
189 >        }
190 >        
191 >        AtomType* atomType;
192 >        if (isDirectional) {
193 >            atomType = new DirectionalAtomType();        
194 >        } else {
195 >            atomType = new AtomType();
196 >        }
197 >        
198          atomType->setName(atomTypeName);
199          atomType->setMass(mass);
200          
201 <        //by default, all of the properties in AtomTypeProperties is set to 0
202 <        //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
203 <        atomType->setLennardJones();
201 >        if (isLJ) {
202 >            atomType->setLennardJones();
203 >        }
204  
205 <        atomType->setIdent(ident);
206 <        atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
207 <        atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
205 >        if (isCharge) {
206 >            atomType->setCharge();
207 >        }
208 >        
209 >        atomType->setIdent(ident);
210 >        
211          atomType->complete();
160 –        //notify a new LJtype atom type is created
161 –        newLJtype(&ident, &sigma, &epsilon, &status);
212  
213 <        //add atom type to AtomTypeContainer
164 <        addAtomType(atomTypeName, atomType);
165 <        ++ident;
213 >        int setLJStatus;
214          
215 <    } else {
216 <        sprintf( painCave.errMsg,
217 <               "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
218 <               "in line %d : %s\n",
219 <               filename.c_str(), lineNo, line);
220 <        painCave.severity = OOPSE_ERROR;
221 <        painCave.isFatal = 1;
222 <        simError();                
215 >        //notify a new LJtype atom type is created
216 >        if (isLJ) {
217 >            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
218 >        }
219 >
220 >        int setChargeStatus;
221 >        if (isCharge) {
222 >            newChargeType(&ident, &charge, &setChargeStatus)
223 >        }
224 >
225 >        if (setLJStatus && setChargeStatus) {
226 >            //add atom type to AtomTypeContainer
227 >            addAtomType(atomTypeName, atomType);
228 >            ++ident;
229 >        } else {
230 >            //error in notifying fortran
231 >            delete atomType;
232 >        }
233      }    
234 +    
235   }
236  
237 +
238 + void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
239 +    StringTokenizer tokenizer(line);
240 +    int nTokens = tokenizer.countTokens();  
241 +    
242 +    //in irectionalAtomTypeSection, a line at least contains 6 tokens
243 +    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
244 +    if (nTokens < 6) {
245 +
246 +    } else {
247 +
248 +
249 +        std::string atomTypeName = tokenizer.nextToken();
250 +        bool isDipole = tokenizer.nextTokenAsBool();
251 +        bool isSticky = tokenizer.nextTokenAsBool();
252 +        double Ixx = tokenizer.nextTokenAsDouble();
253 +        double Iyy = tokenizer.nextTokenAsDouble();
254 +        double Izz = tokenizer.nextTokenAsDouble();
255 +        nTokens -= 6;
256 +
257 +        AtomType* atomType = getAtomType(atomTypeName);
258 +        if (atomType == NULL) {
259 +
260 +        }
261 +        
262 +        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
263 +        if (dAtomType == NULL) {
264 +
265 +
266 +        }
267 +
268 +        if (isDipole) {
269 +            dAtomType->setDipole();
270 +        }
271 +
272 +        if (isSticky) {
273 +            dAtomType->setSticky();
274 +        }
275 +        
276 +        Mat3x3d inertialMat;
277 +        inertialMat(0, 0) = Ixx;
278 +        inertialMat(1, 1) = Ixx;
279 +        inertialMat(2, 2) = Ixx;        
280 +        dAtomType->setI(inertialMat);
281 +
282 +        //read dipole moment
283 +        double dipole;
284 +        if (isDipole) {
285 +            if (nTokens >= 1) {
286 +                dipole = tokenizer.nextTokenAsDouble();
287 +                nTokens -= 1;
288 +            } else {
289 +
290 +            }
291 +        }
292 +
293 +        //read sticky parameters
294 +        double w0;
295 +        double v0;
296 +        double v0p;
297 +        double rl;
298 +        double ru;
299 +        double rlp;
300 +        double rup;        
301 +        if (isSticky) {
302 +            if (nTokens >= 7) {
303 +                w0 = tokenizer.nextTokenAsDouble();
304 +                v0 = tokenizer.nextTokenAsDouble();
305 +                v0p = tokenizer.nextTokenAsDouble();
306 +                rl = tokenizer.nextTokenAsDouble();
307 +                ru = tokenizer.nextTokenAsDouble();
308 +                rlp = tokenizer.nextTokenAsDouble();
309 +                rup = tokenizer.nextTokenAsDouble();    
310 +                nTokens -= 7;
311 +            } else {
312 +
313 +            }
314 +        }
315 +
316 +
317 +        //notify fotran a newDipoleType is created
318 +        int ident = dAtomType->getIdent();
319 +        int setDipoleStatus;
320 +        if (isDipole) {
321 +            newDipoleType(&ident, &dipole, &setDipoleStatus);
322 +        }              
323 +
324 +        //notify fotran a StickyType is created
325 +        int setStickyStatus;
326 +        if (isSticky) {
327 +            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
328 +        }
329 +
330 +
331 +        if (!setDipoleStatus || !setStickyStatus) {
332 +
333 +        }
334 +        
335 +    }
336 + }
337 +
338   void DUFF::parseBondType(const std::string& line, int lineNo){
339  
340      StringTokenizer tokenizer(line);
# Line 256 | Line 416 | void DUFF::parseBondType(const std::string& line, int
416                  for (int i = 0; i < nPairs; ++i) {
417                      power = tokenizer.nextTokenAsInt();
418                      coefficient = tokenizer.nextTokenAsDouble();
419 <
260 <                    if (power < 0) {
261 <
262 <                    } else {
263 <                        pbt->setCoefficient(power, coefficient);
264 <                    }
419 >                    pbt->setCoefficient(power, coefficient);
420                  }
421              }
422              
# Line 374 | Line 529 | void DUFF::parseBendType(const std::string& line, int
529                  for (int i = 0; i < nPairs; ++i) {
530                      power = tokenizer.nextTokenAsInt();
531                      coefficient = tokenizer.nextTokenAsDouble();
532 <
378 <                    if (power < 0) {
379 <
380 <                    } else {
381 <                        pbt->setCoefficient(power, coefficient);
382 <                    }
532 >                    pbt->setCoefficient(power, coefficient);
533                  }
534              }
535              
# Line 463 | Line 613 | void DUFF::parseTorsionType(const std::string& line, i
613                  for (int i = 0; i < nPairs; ++i) {
614                      power = tokenizer.nextTokenAsInt();
615                      coefficient = tokenizer.nextTokenAsDouble();
616 <
467 <                    if (power < 0) {
468 <
469 <                    } else {
470 <                        pbt->setCoefficient(power, coefficient);
471 <                    }
616 >                    pbt->setCoefficient(power, coefficient);
617                  }
618              }
619              
# Line 480 | Line 625 | void DUFF::parseTorsionType(const std::string& line, i
625              } else {
626                  int nSets = nTokens / 3;
627    
483 –                double kchi;
484 –                int n;
485 –                double delta;                
628                  CharmmTorsionType* ctt = new CharmmTorsionType();
629                  
630                  for (int i = 0; i < nSets; ++i) {
631 <                    kchi = tokenizer.nextTokenAsDouble();
632 <                    n = tokenizer.nextTokenAsInt();
633 <                    delta = tokenizer.nextTokenAsDouble();
631 >                    double kchi = tokenizer.nextTokenAsDouble();
632 >                    int n = tokenizer.nextTokenAsInt();
633 >                    double delta = tokenizer.nextTokenAsDouble();
634      
635 <                    if (n < 0 || n > 6) {
494 <
495 <                    } else {
496 <                        ctt->setCharmmTorsionParameter(kchi, n, delta);
497 <                    }
635 >                    ctt->setCharmmTorsionParameter(kchi, n, delta);
636                  }
637              }
638          default:
# Line 507 | Line 645 | void DUFF::parseTorsionType(const std::string& line, i
645   }
646  
647   } //end namespace oopse
510 –

Diff Legend

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