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 1753 by tim, Thu Nov 18 05:42:06 2004 UTC vs.
Revision 1783 by tim, Wed Nov 24 21:33:00 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 + #include "UseTheForce/ForceFieldFactory.hpp"
32 + #include "io/DirectionalAtomTypesSectionParser.hpp"
33 + #include "io/AtomTypesSectionParser.hpp"
34 + #include "io/LennardJonesAtomTypesSectionParser.hpp"
35 + #include "io/ElectrostaticAtomTypesSectionParser.hpp"
36 + #include "io/EAMAtomTypesSectionParser.hpp"
37 + #include "io/StickyAtomTypesSectionParser.hpp"
38 + #include "io/BondTypesSectionParser.hpp"
39 + #include "io/BendTypesSectionParser.hpp"
40 + #include "io/TorsionTypesSectionParser.hpp"
41  
42   namespace oopse {
43  
# Line 33 | Line 46 | ForceField* createDUFF() {
46      return new DUFF();
47   }
48  
49 < //register createDUFF to ForceFieldCreator
50 < ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
49 > //register createDUFF to ForceFieldFactory
50 > bool registerDUFFStatus = ForceFieldFactory::getInstance()->registerForceField("DUFF", createDUFF);
51  
52 + DUFF::DUFF(){
53 +
54 +    //the order of adding section parsers are important
55 +    //DirectionalAtomTypesSectionParser should be added before AtomTypesSectionParser Since
56 +    //These two section parsers will actually create "real" AtomTypes (AtomTypesSectionParser will create
57 +    //AtomType and DirectionalAtomTypesSectionParser will creat DirectionalAtomType which is a subclass
58 +    //of AtomType, therefore it should come first). Other AtomTypes Section Parser will not create the
59 +    //"real" AtomType, they only add and set some attribute of the AtomType. Thus their order are not
60 +    //important. AtomTypesSectionParser should be added before other atom type section parsers.
61 +    //Make sure they are added after DirectionalAtomTypesSectionParser and AtomTypesSectionParser.
62 +    //The order of BondTypesSectionParser, BendTypesSectionParser and TorsionTypesSectionParser are
63 +    //not important.
64 +    spMan_.push_back(new DirectionalAtomTypesSectionParser());
65 +    spMan_.push_back(new AtomTypesSectionParser());
66 +    spMan_.push_back(new LennardJonesAtomTypesSectionParser());
67 +    spMan_.push_back(new ElectrostaticAtomTypesSectionParser());
68 +    spMan_.push_back(new EAMAtomTypesSectionParser());
69 +    spMan_.push_back(new StickyAtomTypesSectionParser());
70 +    spMan_.push_back(new BondTypesSectionParser());
71 +    spMan_.push_back(new BendTypesSectionParser());
72 +    spMan_.push_back(new TorsionTypesSectionParser());
73 +    
74 + }
75 +
76   void DUFF::parse(const std::string& filename) {
77      ifstrstream* ffStream;
78 <    ffStream = openForceFiledFile(filename);
78 >    ffStream = openForceFieldFile(filename);
79 >
80 >    spMan_.parse(*ffStream, *this);
81 >
82 >    ForceField::AtomTypeContainer::MapTypeIterator i;
83 >    AtomType* at;
84 >
85 >    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
86 >        at->makeFortranAtomType();
87 >    }
88 >
89 >    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
90 >        at->complete();
91 >    }
92 >    
93 > }
94 >
95 > /*
96 > ParseState DUFF::getSection(const std::string& section) {
97 >    ParseState result;
98 >    
99 >    switch(section) {
100 >        case "AtomTypes" :
101 >            result = DUFF::AtomTypeSection;
102 >            break;
103 >        case "DirectionalAtomTypes" :
104 >            result = DUFF::DirectionalAtomTypeSection;
105 >            break;
106 >
107 >        case "BondTypes" :
108 >            result = DUFF::BondTypeSection;
109 >            break;
110 >
111 >        case "BendTypes" :
112 >            result = DUFF::BendTypeSection;
113 >            break;
114 >
115 >        case "TorsionTypes" :
116 >            result = DUFF::TorsionTypeSection;
117 >            break;
118 >        default:
119 >            result = DUFF::UnknownSection;
120 >    }
121 >
122 >    return result;
123 > }
124 >
125 > void DUFF::parse(const std::string& filename) {
126 >    ifstrstream* ffStream;
127 >    ffStream = openForceFieldFile(filename);
128      const int bufferSize = 65535;
129 <    char line[bufferSize];
130 <    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)
129 >    std::string line;
130 >    char buffer[bufferSize];
131      int lineNo = 0;
132 <
133 <    while(ffStream.getline(line, bufferSize)){
132 >    int atomIdent = getNAtomType() + 1;  //atom's indent begins from 1 (since only fortran's array begins from 1)
133 >    ParseState currentSection = DUFF::UnknownSection;
134 >    
135 >    while(ffStream.getline(buffer, bufferSize)){
136          ++lineNo;
137  
138 <        //a line begins with '#' or '!' is comment which is skipped
139 <        if (line[0] != '#' ||line[0] != '!') {
140 <            StringTokenizer tokenizer(line);
141 <            
142 <            if (tokenizer.countTokens() >= 4) {
61 <                atomTypeName = tokenizer.nextToken();
62 <                mass = tokenizer.nextTokenAsDouble();
63 <                epsilon = tokenizer.nextTokenAsDouble();
64 <                sigma = tokenizer.nextTokenAsDouble();
138 >        line = trimSpaces(buffer);
139 >        //a line begins with "//" is comment
140 >        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
141 >            continue;
142 >        } else {
143  
144 <                atomType = new AtomType();
145 <                atomType->setName(atomTypeName);
146 <                atomType->setMass(mass);
147 <                
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();
144 >            switch(currentSection) {
145 >                case DUFF::AtomTypeSection :
146 >                    parseAtomType(line, lineNo, atomIdent);
147 >                    break;
148  
149 <                atomType->setIdent(ident);
150 <                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
151 <                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
152 <                atomType->complete();
153 <                //notify a new LJtype atom type is created
154 <                newLJtype(&ident, &sigma, &epsilon, &status);
149 >                case DUFF::DirectionalAtomTypeSection :
150 >                    parseDirectionalAtomType(line, lineNo);
151 >                    break;
152 >                    
153 >                case DUFF::BondTypeSection :
154 >                    parseBondType(line, lineNo);
155 >                    break;
156  
157 <                //add atom type to AtomTypeContainer
158 <                addAtomType(atomTypeName, atomType);
159 <                ++ident;
157 >                case DUFF::BendTypeSection :
158 >                    parseBendType(line, lineNo);
159 >                    break;
160 >
161 >                case DUFF::TorsionTypeSection :
162 >                    parseTorsionType(line, lineNo);
163 >                    break;
164 >                    
165 >                case DUFF::UnknownSection:
166 >                    StringTokenizer tokenizer(line);
167 >                    
168 >                    std::string keyword = tokenizer.nextToken();
169 >                    std::string section = tokenizer.nextToken();
170 >
171 >                    ParseState newSection = getSection(section);
172 >                    if (keyword != "begin" || keyword != "end") {
173 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
174 >                    } else if (keyword == "begin") {
175 >                        if (newSection == DUFF::UnknownSection) {
176 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
177 >                        } else {
178 >                            //enter a new section
179 >                            currentSection = newSection;
180 >                        }
181 >                        
182 >                    } else if (keyword == "end"){
183 >                        if (currentSection == newSection) ) {
184 >                            //leave a section
185 >                            currentSection = DUFF::UnknownSection;
186 >                        } else {                            
187 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
188 >                        }
189 >
190 >                    }
191 >                    break;
192 >                default :
193                  
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();                
194              }
195              
196          }
# Line 98 | Line 199 | void DUFF::parse(const std::string& filename) {
199      delete ffStream;
200   }
201  
202 + void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
203 +    StringTokenizer tokenizer(line);
204 +    int nTokens = tokenizer.countTokens();    
205  
206 < void DUFF::parseAtomType(){
206 >    //in AtomTypeSection, a line at least contains 5 tokens
207 >    //atomTypeName, is Directional, isLJ, isCharge and mass
208 >    if (nTokens < 5)  {
209 >                      
210 >    } else {
211 >
212 >        std::string atomTypeName = tokenizer.nextToken();
213 >        bool isDirectional = tokenizer.nextTokenAsBool();
214 >        bool isLJ = tokenizer.nextTokenAsBool();
215 >        bool isCharge = tokenizer.nextTokenAsBool();        
216 >        double mass = tokenizer.nextTokenAsDouble();
217 >        double epsilon;
218 >        double sigma;
219 >        double charge;
220 >        nTokens -= 5;
221 >
222 >        //parse epsilon and sigma
223 >        if (isLJ) {
224 >            if (nTokens >= 2) {
225 >                epsilon = tokenizer.nextTokenAsDouble();
226 >                sigma = tokenizer.nextTokenAsDouble();
227 >                nTokens -= 2;
228 >            } else {
229 >
230 >            }
231 >        }
232 >
233 >        //parse charge
234 >        if (isCharge) {
235 >            if (nTokens >= 1) {
236 >                charge = tokenizer.nextTokenAsDouble();
237 >                nTokens -= 1;
238 >            } else {
239 >
240 >            }
241 >        }
242 >        
243 >        AtomType* atomType;
244 >        if (isDirectional) {
245 >            atomType = new DirectionalAtomType();        
246 >        } else {
247 >            atomType = new AtomType();
248 >        }
249 >        
250 >        atomType->setName(atomTypeName);
251 >        atomType->setMass(mass);
252 >        
253 >        if (isLJ) {
254 >            atomType->setLennardJones();
255 >        }
256 >
257 >        if (isCharge) {
258 >            atomType->setCharge();
259 >        }
260 >        
261 >        atomType->setIdent(ident);
262 >        
263 >        atomType->complete();
264 >
265 >        int setLJStatus;
266 >        
267 >        //notify a new LJtype atom type is created
268 >        if (isLJ) {
269 >            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
270 >        }
271 >
272 >        int setChargeStatus;
273 >        if (isCharge) {
274 >            newChargeType(&ident, &charge, &setChargeStatus)
275 >        }
276 >
277 >        if (setLJStatus && setChargeStatus) {
278 >            //add atom type to AtomTypeContainer
279 >            addAtomType(atomTypeName, atomType);
280 >            ++ident;
281 >        } else {
282 >            //error in notifying fortran
283 >            delete atomType;
284 >        }
285 >    }    
286 >    
287   }
288  
289 < void DUFF::parseBondType(){
289 >
290 > void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
291 >    StringTokenizer tokenizer(line);
292 >    int nTokens = tokenizer.countTokens();  
293 >    
294 >    //in DirectionalAtomTypeSection, a line at least contains 6 tokens
295 >    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
296 >    if (nTokens < 6) {
297 >        std::cerr << "Not enought tokens" << std::endl;
298 >    } else {
299 >
300 >
301 >        std::string atomTypeName = tokenizer.nextToken();
302 >        bool isDipole = tokenizer.nextTokenAsBool();
303 >        bool isSticky = tokenizer.nextTokenAsBool();
304 >        double Ixx = tokenizer.nextTokenAsDouble();
305 >        double Iyy = tokenizer.nextTokenAsDouble();
306 >        double Izz = tokenizer.nextTokenAsDouble();
307 >        nTokens -= 6;
308 >
309 >        AtomType* atomType = getAtomType(atomTypeName);
310 >        if (atomType == NULL) {
311 >
312 >        }
313 >        
314 >        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
315 >        if (dAtomType == NULL) {
316 >
317 >
318 >        }
319 >
320 >        if (isDipole) {
321 >            dAtomType->setDipole();
322 >        }
323 >
324 >        if (isSticky) {
325 >            dAtomType->setSticky();
326 >        }
327 >        
328 >        Mat3x3d inertialMat;
329 >        inertialMat(0, 0) = Ixx;
330 >        inertialMat(1, 1) = Iyy;
331 >        inertialMat(2, 2) = Izz;        
332 >        dAtomType->setI(inertialMat);
333 >
334 >        //read dipole moment
335 >        double dipole;
336 >        if (isDipole) {
337 >            if (nTokens >= 1) {
338 >                dipole = tokenizer.nextTokenAsDouble();
339 >                nTokens -= 1;
340 >            } else {
341 >
342 >            }
343 >        }
344 >
345 >        //read sticky parameters
346 >        double w0;
347 >        double v0;
348 >        double v0p;
349 >        double rl;
350 >        double ru;
351 >        double rlp;
352 >        double rup;        
353 >        if (isSticky) {
354 >            if (nTokens >= 7) {
355 >                w0 = tokenizer.nextTokenAsDouble();
356 >                v0 = tokenizer.nextTokenAsDouble();
357 >                v0p = tokenizer.nextTokenAsDouble();
358 >                rl = tokenizer.nextTokenAsDouble();
359 >                ru = tokenizer.nextTokenAsDouble();
360 >                rlp = tokenizer.nextTokenAsDouble();
361 >                rup = tokenizer.nextTokenAsDouble();    
362 >                nTokens -= 7;
363 >            } else {
364 >
365 >            }
366 >        }
367 >
368 >
369 >        //notify fotran a newDipoleType is created
370 >        int ident = dAtomType->getIdent();
371 >        int setDipoleStatus;
372 >        if (isDipole) {
373 >            newDipoleType(&ident, &dipole, &setDipoleStatus);
374 >        }              
375 >
376 >        //notify fotran a StickyType is created
377 >        int setStickyStatus;
378 >        if (isSticky) {
379 >            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
380 >        }
381 >
382 >
383 >        if (!setDipoleStatus || !setStickyStatus) {
384 >
385 >        }
386 >        
387 >    }
388   }
389  
390 < void DUFF::parseBendType(){
390 > void DUFF::parseBondType(const std::string& line, int lineNo){
391 >
392 >    StringTokenizer tokenizer(line);
393 >    std::string at1;
394 >    std::string at2;
395 >    std::string bt;
396 >    BondType* bondType = NULL;
397 >    double b0;
398 >    
399 >    int nTokens = tokenizer.countTokens();
400 >
401 >    if (nTokens < 4) {
402 >
403 >        return;
404 >    }
405 >    
406 >    at1 = tokenizer.nextToken();
407 >    at2 = tokenizer.nextToken();
408 >    bt = tokenizer.nextToken();
409 >    b0 = tokenizer.nextTokenAsDouble();
410 >    nTokens -= 4;
411 >
412 >    //switch is a maintain nightmare
413 >    switch(bt) {
414 >        case "Fixed" :
415 >            bondType = new FixedBondType();
416 >            break;
417 >            
418 >        case "Harmonic" :
419 >            if (nTokens < 1) {
420 >
421 >            } else {
422 >
423 >                double kb = tokenizer.nextTokenAsDouble();
424 >                bondType = new HarmonicBondType(b0, kb);
425 >            }
426 >
427 >            break;
428 >
429 >        case "Cubic" :
430 >            if (nTokens < 4) {
431 >
432 >            } else {
433 >
434 >                double k3 = tokenizer.nextTokenAsDouble();
435 >                double k2 = tokenizer.nextTokenAsDouble();
436 >                double k1 = tokenizer.nextTokenAsDouble();
437 >                double k0 = tokenizer.nextTokenAsDouble();
438 >                
439 >                bondType = new CubicBondType(b0, k3, k2, k1, k0);
440 >            }
441 >            break;
442 >            
443 >        case "Quartic" :
444 >            if (nTokens < 5) {
445 >
446 >            } else {
447 >
448 >                b0 = tokenizer.nextTokenAsDouble();
449 >                double k4 = tokenizer.nextTokenAsDouble();
450 >                double k3 = tokenizer.nextTokenAsDouble();
451 >                double k2 = tokenizer.nextTokenAsDouble();
452 >                double k1 = tokenizer.nextTokenAsDouble();
453 >                double k0 = tokenizer.nextTokenAsDouble();
454 >                
455 >                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
456 >            }
457 >            break;
458 >
459 >        case "Polynomial" :
460 >            if (nTokens < 2 || nTokens % 2 != 0) {
461 >
462 >            } else {
463 >                int nPairs = nTokens / 2;
464 >                int power;
465 >                double coefficient;
466 >                PolynomialBondType pbt = new PolynomialBondType();
467 >                
468 >                for (int i = 0; i < nPairs; ++i) {
469 >                    power = tokenizer.nextTokenAsInt();
470 >                    coefficient = tokenizer.nextTokenAsDouble();
471 >                    pbt->setCoefficient(power, coefficient);
472 >                }
473 >            }
474 >            
475 >            break;
476 >
477 >        default:
478 >            
479 >    }
480 >
481 >    if (bondType != NULL) {
482 >        addBondType(at1, at2, bondType);
483 >    }
484   }
485  
486 < void DUFF::parseTorsionType(){
486 > void DUFF::parseBendType(const std::string& line, int lineNo){
487 >    StringTokenizer tokenizer(line);
488 >    std::string at1;
489 >    std::string at2;
490 >    std::string at3;
491 >    std::string bt;
492 >    double theta0;
493 >    BendType* bendType = NULL;
494  
495 +    int nTokens = tokenizer.countTokens();
496 +
497 +    if (nTokens < 5) {
498 +
499 +        return;
500 +    }
501 +    
502 +    at1 = tokenizer.nextToken();
503 +    at2 = tokenizer.nextToken();
504 +    at3 = tokenizer.nextToken();
505 +    bt = tokenizer.nextToken();
506 +    theta0 = tokenizer.nextTokenAsDouble();
507 +    nTokens -= 5;
508 +
509 +    //switch is a maintain nightmare
510 +    switch(bt) {
511 +            
512 +        case "Harmonic" :
513 +            
514 +            if (nTokens < 1) {
515 +
516 +            } else {
517 +
518 +                double ktheta = tokenizer.nextTokenAsDouble();
519 +                bendType = new HarmonicBendType(theta0, ktheta);
520 +            }
521 +            break;
522 +        case "GhostBend" :
523 +            if (nTokens < 1) {
524 +
525 +            } else {
526 +                double ktheta = tokenizer.nextTokenAsDouble();
527 +                bendType = new HarmonicBendType(theta0, ktheta);                
528 +            }
529 +            break;            
530 +
531 +        case "UreyBradley" :
532 +            if (nTokens < 3) {
533 +
534 +            } else {
535 +                double ktheta = tokenizer.nextTokenAsDouble();
536 +                double s0 =  tokenizer.nextTokenAsDouble();
537 +                double kub = tokenizer.nextTokenAsDouble();
538 +                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
539 +            }
540 +            break;
541 +            
542 +        case "Cubic" :
543 +            if (nTokens < 4) {
544 +
545 +            } else {
546 +
547 +                double k3 = tokenizer.nextTokenAsDouble();
548 +                double k2 = tokenizer.nextTokenAsDouble();
549 +                double k1 = tokenizer.nextTokenAsDouble();
550 +                double k0 = tokenizer.nextTokenAsDouble();
551 +                
552 +                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
553 +            }
554 +            break;
555 +            
556 +        case "Quartic" :
557 +            if (nTokens < 5) {
558 +
559 +            } else {
560 +
561 +                theta0 = tokenizer.nextTokenAsDouble();
562 +                double k4 = tokenizer.nextTokenAsDouble();
563 +                double k3 = tokenizer.nextTokenAsDouble();
564 +                double k2 = tokenizer.nextTokenAsDouble();
565 +                double k1 = tokenizer.nextTokenAsDouble();
566 +                double k0 = tokenizer.nextTokenAsDouble();
567 +                
568 +                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
569 +            }
570 +            break;
571 +
572 +        case "Polynomial" :
573 +            if (nTokens < 2 || nTokens % 2 != 0) {
574 +
575 +            } else {
576 +                int nPairs = nTokens / 2;
577 +                int power;
578 +                double coefficient;
579 +                PolynomialBendType* pbt = new PolynomialBendType();
580 +                
581 +                for (int i = 0; i < nPairs; ++i) {
582 +                    power = tokenizer.nextTokenAsInt();
583 +                    coefficient = tokenizer.nextTokenAsDouble();
584 +                    pbt->setCoefficient(power, coefficient);
585 +                }
586 +            }
587 +            
588 +            break;
589 +
590 +        default:
591 +            
592 +    }
593 +
594 +    if (bendType != NULL) {
595 +        addBendType(at1, at2, at3, bendType);
596 +    }
597 +    
598   }
599  
600 < } //end namespace oopse
600 > void DUFF::parseTorsionType(const std::string& line, int lineNo){
601 >    StringTokenizer tokenizer(line);
602 >    std::string at1;
603 >    std::string at2;
604 >    std::string at3;
605 >    std::string at4;
606 >    std::string tt;
607 >    TorsionType* torsionType = NULL;
608  
609 +    int nTokens = tokenizer.countTokens();
610 +
611 +    if (nTokens < 5) {
612 +
613 +        return;
614 +    }
615 +    
616 +    at1 = tokenizer.nextToken();
617 +    at2 = tokenizer.nextToken();
618 +    at3 = tokenizer.nextToken();
619 +    at4 = tokenizer.nextToken();
620 +    tt = tokenizer.nextToken();
621 +
622 +    nTokens -= 5;
623 +
624 +    switch(tt) {
625 +            
626 +        case "Cubic" :
627 +            if (nTokens < 4) {
628 +
629 +            } else {
630 +
631 +                double k3 = tokenizer.nextTokenAsDouble();
632 +                double k2 = tokenizer.nextTokenAsDouble();
633 +                double k1 = tokenizer.nextTokenAsDouble();
634 +                double k0 = tokenizer.nextTokenAsDouble();
635 +                
636 +                bendType = new CubicTorsionType(k3, k2, k1, k0);
637 +            }
638 +            break;
639 +            
640 +        case "Quartic" :
641 +            if (nTokens < 5) {
642 +
643 +            } else {
644 +
645 +                theta0 = tokenizer.nextTokenAsDouble();
646 +                double k4 = tokenizer.nextTokenAsDouble();
647 +                double k3 = tokenizer.nextTokenAsDouble();
648 +                double k2 = tokenizer.nextTokenAsDouble();
649 +                double k1 = tokenizer.nextTokenAsDouble();
650 +                double k0 = tokenizer.nextTokenAsDouble();
651 +                
652 +                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
653 +            }
654 +            break;
655 +
656 +        case "Polynomial" :
657 +            if (nTokens < 2 || nTokens % 2 != 0) {
658 +
659 +            } else {
660 +                int nPairs = nTokens / 2;
661 +                int power;
662 +                double coefficient;
663 +                PolynomialTorsionType* pbt = new PolynomialTorsionType();
664 +                
665 +                for (int i = 0; i < nPairs; ++i) {
666 +                    power = tokenizer.nextTokenAsInt();
667 +                    coefficient = tokenizer.nextTokenAsDouble();
668 +                    pbt->setCoefficient(power, coefficient);
669 +                }
670 +            }
671 +            
672 +            break;
673 +        case "Charmm" :
674 +            
675 +            if (nTokens < 3 || nTokens % 3 != 0) {
676 +
677 +            } else {
678 +                int nSets = nTokens / 3;
679 +  
680 +                CharmmTorsionType* ctt = new CharmmTorsionType();
681 +                
682 +                for (int i = 0; i < nSets; ++i) {
683 +                    double kchi = tokenizer.nextTokenAsDouble();
684 +                    int n = tokenizer.nextTokenAsInt();
685 +                    double delta = tokenizer.nextTokenAsDouble();
686 +    
687 +                    ctt->setCharmmTorsionParameter(kchi, n, delta);
688 +                }
689 +            }
690 +        default:
691 +            
692 +    }
693 +
694 +    if (torsionType != NULL) {
695 +        addTorsionType(at1, at2, at3, at4, torsionType);
696 +    }
697 + }
698 + */
699 + } //end namespace oopse

Diff Legend

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