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 1770 by tim, Tue Nov 23 17:53:43 2004 UTC vs.
Revision 1840 by tim, Fri Dec 3 00:26:07 2004 UTC

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

Diff Legend

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