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 1756 by tim, Thu Nov 18 05:42:06 2004 UTC vs.
Revision 1757 by tim, Fri Nov 19 00:11:33 2004 UTC

# Line 36 | Line 36 | void DUFF::parse(const std::string& filename) {
36   //register createDUFF to ForceFieldCreator
37   ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
38  
39 +
40 + ParseState DUFF::getSection(const std::string& section) {
41 +    ParseState result;
42 +    
43 +    switch(section) {
44 +        case "AtomTypes" :
45 +            result = DUFF::AtomTypeSection;
46 +            break;
47 +
48 +        case "BondTypes" :
49 +            result = DUFF::BondTypeSection;
50 +            break;
51 +
52 +        case "BendTypes" :
53 +            result = DUFF::BendTypeSection;
54 +            break;
55 +
56 +        case "TorsionTypes" :
57 +            result = DUFF::TorsionTypeSection;
58 +            break;
59 +        default:
60 +            result = DUFF::UnknownSection;
61 +    }
62 +
63 +    return result;
64 + }
65 +
66   void DUFF::parse(const std::string& filename) {
67      ifstrstream* ffStream;
68      ffStream = openForceFiledFile(filename);
69      const int bufferSize = 65535;
70 <    char line[bufferSize];
71 <    AtomType* atomType;
45 <    std::string atomTypeName;
46 <    double mass;
47 <    double epsilon;
48 <    double sigma;
49 <    int status;
50 <    int ident = 1; //ident begins from 1 (fortran's index begins from 1)
70 >    std::string line;
71 >    char buffer[bufferSize];
72      int lineNo = 0;
73 <
74 <    while(ffStream.getline(line, bufferSize)){
73 >    ParseState currentSection = DUFF::UnknownSection;
74 >    
75 >    while(ffStream.getline(buffer, bufferSize)){
76          ++lineNo;
77  
78 <        //a line begins with '#' or '!' is comment which is skipped
79 <        if (line[0] != '#' ||line[0] != '!') {
80 <            StringTokenizer tokenizer(line);
81 <            
82 <            if (tokenizer.countTokens() >= 4) {
61 <                atomTypeName = tokenizer.nextToken();
62 <                mass = tokenizer.nextTokenAsDouble();
63 <                epsilon = tokenizer.nextTokenAsDouble();
64 <                sigma = tokenizer.nextTokenAsDouble();
78 >        line = trimSpaces(buffer);
79 >        //a line begins with "//" is comment
80 >        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
81 >            continue;
82 >        } else {
83  
84 <                atomType = new AtomType();
85 <                atomType->setName(atomTypeName);
86 <                atomType->setMass(mass);
87 <                
70 <                //by default, all of the properties in AtomTypeProperties is set to 0
71 <                //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
72 <                atomType->setLennardJones();
84 >            switch(currentSection) {
85 >                case DUFF::AtomTypeSection :
86 >                    parseAtomType(line, lineNo);
87 >                    break;
88  
89 <                atomType->setIdent(ident);
90 <                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
91 <                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
77 <                atomType->complete();
78 <                //notify a new LJtype atom type is created
79 <                newLJtype(&ident, &sigma, &epsilon, &status);
89 >                case DUFF::BondTypeSection :
90 >                    parseBondType(line, lineNo);
91 >                    break;
92  
93 <                //add atom type to AtomTypeContainer
94 <                addAtomType(atomTypeName, atomType);
95 <                ++ident;
93 >                case DUFF::BendTypeSection :
94 >                    parseBendType(line, lineNo);
95 >                    break;
96 >
97 >                case DUFF::TorsionTypeSection :
98 >                    parseTorsionType(line, lineNo);
99 >                    break;
100 >                    
101 >                case DUFF::UnknownSection:
102 >                    StringTokenizer tokenizer(line);
103 >                    
104 >                    std::string keyword = tokenizer.nextToken();
105 >                    std::string section = tokenizer.nextToken();
106 >
107 >                    ParseState newSection = getSection(section);
108 >                    if (keyword != "begin" || keyword != "end") {
109 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
110 >                    } else if (keyword == "begin") {
111 >                        if (newSection == DUFF::UnknownSection) {
112 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
113 >                        } else {
114 >                            //enter a new section
115 >                            currentSection = newSection;
116 >                        }
117 >                        
118 >                    } else if (keyword == "end"){
119 >                        if (currentSection == newSection) ) {
120 >                            //leave a section
121 >                            currentSection = DUFF::UnknownSection;
122 >                        } else {                            
123 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
124 >                        }
125 >
126 >                    }
127 >                    break;
128 >                default :
129                  
85 –            } else {
86 –                sprintf( painCave.errMsg,
87 –                       "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
88 –                       "in line %d : %s\n",
89 –                       filename.c_str(), lineNo, line);
90 –                painCave.severity = OOPSE_ERROR;
91 –                painCave.isFatal = 1;
92 –                simError();                
130              }
131              
132          }
# Line 99 | Line 136 | void DUFF::parseAtomType(){
136   }
137  
138  
139 < void DUFF::parseAtomType(){
139 > void DUFF::parseAtomType(const std::string& line, int lineNo){
140 >    StringTokenizer tokenizer(line);
141 >    
142 >    if (tokenizer.countTokens() >= 4) {
143 >        atomTypeName = tokenizer.nextToken();
144 >        mass = tokenizer.nextTokenAsFloat();
145 >        epsilon = tokenizer.nextTokenAsFloat();
146 >        sigma = tokenizer.nextTokenAsFloat();
147 >
148 >        atomType = new AtomType();
149 >        atomType->setName(atomTypeName);
150 >        atomType->setMass(mass);
151 >        
152 >        //by default, all of the properties in AtomTypeProperties is set to 0
153 >        //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
154 >        atomType->setLennardJones();
155 >
156 >        atomType->setIdent(ident);
157 >        atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
158 >        atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
159 >        atomType->complete();
160 >        //notify a new LJtype atom type is created
161 >        newLJtype(&ident, &sigma, &epsilon, &status);
162 >
163 >        //add atom type to AtomTypeContainer
164 >        addAtomType(atomTypeName, atomType);
165 >        ++ident;
166 >        
167 >    } else {
168 >        sprintf( painCave.errMsg,
169 >               "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
170 >               "in line %d : %s\n",
171 >               filename.c_str(), lineNo, line);
172 >        painCave.severity = OOPSE_ERROR;
173 >        painCave.isFatal = 1;
174 >        simError();                
175 >    }    
176   }
177  
178 < void DUFF::parseBondType(){
178 > void DUFF::parseBondType(const std::string& line, int lineNo){
179 >
180 >    StringTokenizer tokenizer(line);
181 >    std::string at1;
182 >    std::string at2;
183 >    std::string bt;
184 >    BondType* bondType = NULL;
185 >    double b0;
186 >    
187 >    int nTokens = tokenizer.countTokens();
188 >
189 >    if (nTokens < 4) {
190 >
191 >        return;
192 >    }
193 >    
194 >    at1 = tokenizer.nextToken();
195 >    at2 = tokenizer.nextToken();
196 >    bt = tokenizer.nextToken();
197 >    b0 = tokenizer.nextTokenAsDouble();
198 >    nTokens -= 4;
199 >
200 >    //switch is a maintain nightmare
201 >    switch(bt) {
202 >        case "FixedBondType" :
203 >            bondType = new FixedBondType();
204 >            break;
205 >            
206 >        case "HarmonicBondType" :
207 >            if (nTokens < 1) {
208 >
209 >            } else {
210 >
211 >                double kb = tokenizer.nextTokenAsDouble();
212 >                bondType = new HarmonicBondType(b0, kb);
213 >            }
214 >
215 >            break;
216 >
217 >        case "CubicBondType" :
218 >            if (nTokens < 4) {
219 >
220 >            } else {
221 >
222 >                double k3 = tokenizer.nextTokenAsDouble();
223 >                double k2 = tokenizer.nextTokenAsDouble();
224 >                double k1 = tokenizer.nextTokenAsDouble();
225 >                double k0 = tokenizer.nextTokenAsDouble();
226 >                
227 >                bondType = new CubicBondType(b0, k3, k2, k1, k0);
228 >            }
229 >            break;
230 >            
231 >        case "QuadraticBondType" :
232 >            if (nTokens < 5) {
233 >
234 >            } else {
235 >
236 >                b0 = tokenizer.nextTokenAsDouble();
237 >                double k4 = tokenizer.nextTokenAsDouble();
238 >                double k3 = tokenizer.nextTokenAsDouble();
239 >                double k2 = tokenizer.nextTokenAsDouble();
240 >                double k1 = tokenizer.nextTokenAsDouble();
241 >                double k0 = tokenizer.nextTokenAsDouble();
242 >                
243 >                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
244 >            }
245 >            break;
246 >
247 >        case "PolynomialBondType " :
248 >            if (nTokens < 2 || nTokens % 2 != 0) {
249 >
250 >            } else {
251 >                int nPairs = nTokens / 2;
252 >                int power;
253 >                double coefficient;
254 >                PolynomialBondType pbt = new PolynomialBondType();
255 >                
256 >                for (int i = 0; i < nPairs; ++i) {
257 >                    power = tokenizer.nextTokenAsInt();
258 >                    coefficient = tokenizer.nextTokenAsDouble();
259 >
260 >                    if (power < 0) {
261 >
262 >                    } else {
263 >                        pbt->setCoefficient(power, coefficient);
264 >                    }
265 >                }
266 >            }
267 >            
268 >            break;
269 >
270 >        default:
271 >            
272 >    }
273 >
274 >    if (bondType != NULL) {
275 >        addBondType(at1, at2, bondType);
276 >    }
277   }
278  
279 < void DUFF::parseBendType(){
279 > void DUFF::parseBendType(const std::string& line, int lineNo){
280 >    StringTokenizer tokenizer(line);
281 >    std::string at1;
282 >    std::string at2;
283 >    std::string at3;
284 >    std::string bt;
285 >    double theta0;
286 >    BendType* bendType = NULL;
287 >
288 >    int nTokens = tokenizer.countTokens();
289 >
290 >    if (nTokens < 5) {
291 >
292 >        return;
293 >    }
294 >    
295 >    at1 = tokenizer.nextToken();
296 >    at2 = tokenizer.nextToken();
297 >    at3 = tokenizer.nextToken();
298 >    bt = tokenizer.nextToken();
299 >    theta0 = tokenizer.nextTokenAsDouble();
300 >    nTokens -= 5;
301 >
302 >    //switch is a maintain nightmare
303 >    switch(bt) {
304 >            
305 >        case "HarmonicBendType" :
306 >            
307 >            if (nTokens < 1) {
308 >
309 >            } else {
310 >
311 >                double ktheta = tokenizer.nextTokenAsDouble();
312 >                bendType = new HarmonicBendType(theta0, ktheta);
313 >            }
314 >            break;
315 >        case "GhostBendType" :
316 >            if (nTokens < 1) {
317 >
318 >            } else {
319 >                double ktheta = tokenizer.nextTokenAsDouble();
320 >                bendType = new HarmonicBendType(theta0, ktheta);                
321 >            }
322 >            break;            
323 >
324 >        case "UreyBradleyBendType" :
325 >            if (nTokens < 3) {
326 >
327 >            } else {
328 >                double ktheta = tokenizer.nextTokenAsDouble();
329 >                double s0 =  tokenizer.nextTokenAsDouble();
330 >                double kub = tokenizer.nextTokenAsDouble();
331 >                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
332 >            }
333 >            break;
334 >            
335 >        case "CubicBendType" :
336 >            if (nTokens < 4) {
337 >
338 >            } else {
339 >
340 >                double k3 = tokenizer.nextTokenAsDouble();
341 >                double k2 = tokenizer.nextTokenAsDouble();
342 >                double k1 = tokenizer.nextTokenAsDouble();
343 >                double k0 = tokenizer.nextTokenAsDouble();
344 >                
345 >                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
346 >            }
347 >            break;
348 >            
349 >        case "QuadraticBendType" :
350 >            if (nTokens < 5) {
351 >
352 >            } else {
353 >
354 >                theta0 = tokenizer.nextTokenAsDouble();
355 >                double k4 = tokenizer.nextTokenAsDouble();
356 >                double k3 = tokenizer.nextTokenAsDouble();
357 >                double k2 = tokenizer.nextTokenAsDouble();
358 >                double k1 = tokenizer.nextTokenAsDouble();
359 >                double k0 = tokenizer.nextTokenAsDouble();
360 >                
361 >                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
362 >            }
363 >            break;
364 >
365 >        case "PolynomialBendType " :
366 >            if (nTokens < 2 || nTokens % 2 != 0) {
367 >
368 >            } else {
369 >                int nPairs = nTokens / 2;
370 >                int power;
371 >                double coefficient;
372 >                PolynomialBendType* pbt = new PolynomialBendType();
373 >                
374 >                for (int i = 0; i < nPairs; ++i) {
375 >                    power = tokenizer.nextTokenAsInt();
376 >                    coefficient = tokenizer.nextTokenAsDouble();
377 >
378 >                    if (power < 0) {
379 >
380 >                    } else {
381 >                        pbt->setCoefficient(power, coefficient);
382 >                    }
383 >                }
384 >            }
385 >            
386 >            break;
387 >
388 >        default:
389 >            
390 >    }
391 >
392 >    if (bendType != NULL) {
393 >        addBendType(at1, at2, at3, bendType);
394 >    }
395 >    
396   }
397  
398 < void DUFF::parseTorsionType(){
398 > void DUFF::parseTorsionType(const std::string& line, int lineNo){
399 >    StringTokenizer tokenizer(line);
400 >    std::string at1;
401 >    std::string at2;
402 >    std::string at3;
403 >    std::string at4;
404 >    std::string tt;
405 >    TorsionType* bendType = NULL;
406  
407 +    int nTokens = tokenizer.countTokens();
408 +
409 +    if (nTokens < 5) {
410 +
411 +        return;
412 +    }
413 +    
414 +    at1 = tokenizer.nextToken();
415 +    at2 = tokenizer.nextToken();
416 +    at3 = tokenizer.nextToken();
417 +    at4 = tokenizer.nextToken();
418 +    tt = tokenizer.nextToken();
419 +
420 +    nTokens -= 5;
421 +
422 +    switch(tt) {
423 +            
424 +        case "CubicTorsionType" :
425 +            if (nTokens < 4) {
426 +
427 +            } else {
428 +
429 +                double k3 = tokenizer.nextTokenAsDouble();
430 +                double k2 = tokenizer.nextTokenAsDouble();
431 +                double k1 = tokenizer.nextTokenAsDouble();
432 +                double k0 = tokenizer.nextTokenAsDouble();
433 +                
434 +                bendType = new CubicTorsionType(k3, k2, k1, k0);
435 +            }
436 +            break;
437 +            
438 +        case "QuadraticTorsionType" :
439 +            if (nTokens < 5) {
440 +
441 +            } else {
442 +
443 +                theta0 = tokenizer.nextTokenAsDouble();
444 +                double k4 = tokenizer.nextTokenAsDouble();
445 +                double k3 = tokenizer.nextTokenAsDouble();
446 +                double k2 = tokenizer.nextTokenAsDouble();
447 +                double k1 = tokenizer.nextTokenAsDouble();
448 +                double k0 = tokenizer.nextTokenAsDouble();
449 +                
450 +                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
451 +            }
452 +            break;
453 +
454 +        case "PolynomialTorsionType " :
455 +            if (nTokens < 2 || nTokens % 2 != 0) {
456 +
457 +            } else {
458 +                int nPairs = nTokens / 2;
459 +                int power;
460 +                double coefficient;
461 +                PolynomialTorsionType* pbt = new PolynomialTorsionType();
462 +                
463 +                for (int i = 0; i < nPairs; ++i) {
464 +                    power = tokenizer.nextTokenAsInt();
465 +                    coefficient = tokenizer.nextTokenAsDouble();
466 +
467 +                    if (power < 0) {
468 +
469 +                    } else {
470 +                        pbt->setCoefficient(power, coefficient);
471 +                    }
472 +                }
473 +            }
474 +            
475 +            break;
476 +        case "CharmmTorsionType" :
477 +            
478 +            if (nTokens < 3 || nTokens % 3 != 0) {
479 +
480 +            } else {
481 +                int nSets = nTokens / 3;
482 +  
483 +                double kchi;
484 +                int n;
485 +                double delta;                
486 +                CharmmTorsionType* ctt = new CharmmTorsionType();
487 +                
488 +                for (int i = 0; i < nSets; ++i) {
489 +                    kchi = tokenizer.nextTokenAsDouble();
490 +                    n = tokenizer.nextTokenAsInt();
491 +                    delta = tokenizer.nextTokenAsDouble();
492 +    
493 +                    if (n < 0 || n > 6) {
494 +
495 +                    } else {
496 +                        ctt->setCharmmTorsionParameter(kchi, n, delta);
497 +                    }
498 +                }
499 +            }
500 +        default:
501 +            
502 +    }
503 +
504 +    if (bendType != NULL) {
505 +        addTorsionType(at1, at2, at3, bendType);
506 +    }
507   }
508  
509   } //end namespace oopse

Diff Legend

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