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 1767 by tim, Mon Nov 22 23:04:16 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  
32   namespace oopse {
33  
# Line 33 | Line 36 | ForceField* createDUFF() {
36      return new DUFF();
37   }
38  
39 < //register createDUFF to ForceFieldCreator
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 +    
46 +    switch(section) {
47 +        case "AtomTypes" :
48 +            result = DUFF::AtomTypeSection;
49 +            break;
50 +        case "DirectionalAtomTypes" :
51 +            result = DUFF::DirectionalAtomTypeSection;
52 +            break;
53 +
54 +        case "BondTypes" :
55 +            result = DUFF::BondTypeSection;
56 +            break;
57 +
58 +        case "BendTypes" :
59 +            result = DUFF::BendTypeSection;
60 +            break;
61 +
62 +        case "TorsionTypes" :
63 +            result = DUFF::TorsionTypeSection;
64 +            break;
65 +        default:
66 +            result = DUFF::UnknownSection;
67 +    }
68 +
69 +    return result;
70 + }
71 +
72   void DUFF::parse(const std::string& filename) {
73      ifstrstream* ffStream;
74      ffStream = openForceFiledFile(filename);
75      const int bufferSize = 65535;
76 <    char line[bufferSize];
77 <    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)
76 >    std::string line;
77 >    char buffer[bufferSize];
78      int lineNo = 0;
79 <
80 <    while(ffStream.getline(line, bufferSize)){
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;
84  
85 <        //a line begins with '#' or '!' is comment which is skipped
86 <        if (line[0] != '#' ||line[0] != '!') {
87 <            StringTokenizer tokenizer(line);
88 <            
89 <            if (tokenizer.countTokens() >= 4) {
61 <                atomTypeName = tokenizer.nextToken();
62 <                mass = tokenizer.nextTokenAsDouble();
63 <                epsilon = tokenizer.nextTokenAsDouble();
64 <                sigma = tokenizer.nextTokenAsDouble();
85 >        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 {
90  
91 <                atomType = new AtomType();
92 <                atomType->setName(atomTypeName);
93 <                atomType->setMass(mass);
94 <                
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();
91 >            switch(currentSection) {
92 >                case DUFF::AtomTypeSection :
93 >                    parseAtomType(line, lineNo, atomIdent);
94 >                    break;
95  
96 <                atomType->setIdent(ident);
97 <                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
98 <                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
99 <                atomType->complete();
100 <                //notify a new LJtype atom type is created
101 <                newLJtype(&ident, &sigma, &epsilon, &status);
96 >                case DUFF::DirectionalAtomTypeSection :
97 >                    parseDirectionalAtomType(line, lineNo);
98 >                    break;
99 >                    
100 >                case DUFF::BondTypeSection :
101 >                    parseBondType(line, lineNo);
102 >                    break;
103  
104 <                //add atom type to AtomTypeContainer
105 <                addAtomType(atomTypeName, atomType);
106 <                ++ident;
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                  
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();                
141              }
142              
143          }
# Line 98 | Line 146 | void DUFF::parse(const std::string& filename) {
146      delete ffStream;
147   }
148  
149 + void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
150 +    StringTokenizer tokenizer(line);
151 +    int nTokens = tokenizer.countTokens();    
152  
153 < void DUFF::parseAtomType(){
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 < void DUFF::parseBondType(){
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 >        
334 >    }
335   }
336  
337 < void DUFF::parseBendType(){
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;
345 >    
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 "Fixed" :
362 >            bondType = new FixedBondType();
363 >            break;
364 >            
365 >        case "Harmonic" :
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 "Cubic" :
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 "Quartic" :
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 "Polynomial" :
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 >    }
431   }
432  
433 < void DUFF::parseTorsionType(){
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 "Harmonic" :
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 "GhostBend" :
470 +            if (nTokens < 1) {
471 +
472 +            } else {
473 +                double ktheta = tokenizer.nextTokenAsDouble();
474 +                bendType = new HarmonicBendType(theta0, ktheta);                
475 +            }
476 +            break;            
477 +
478 +        case "UreyBradley" :
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 "Cubic" :
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 "Quartic" :
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 "Polynomial" :
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 < } //end namespace oopse
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* torsionType = 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 "Cubic" :
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 "Quartic" :
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 "Polynomial" :
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 "Charmm" :
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 (torsionType != NULL) {
642 +        addTorsionType(at1, at2, at3, at4, torsionType);
643 +    }
644 + }
645 +
646 + } //end namespace oopse

Diff Legend

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