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

Comparing branches/new_design/OOPSE-4/src/UseTheForce/DUFF.cpp (file contents):
Revision 1836 by tim, Tue Nov 30 22:43:51 2004 UTC vs.
Revision 1837 by tim, Thu Dec 2 22:15:31 2004 UTC

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

Diff Legend

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