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

Diff Legend

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