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 1760 by tim, Fri Nov 19 20:23:26 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 99 | Line 147 | void DUFF::parseAtomType(){
147   }
148  
149  
150 < void DUFF::parseAtomType(){
150 > void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
151 >    StringTokenizer tokenizer(line);
152 >    int nTokens = tokenizer.countTokens();    
153 >
154 >    //in AtomTypeSection, a line at least contains 5 tokens
155 >    //atomTypeName, is Directional, isLJ, isCharge and mass
156 >    if (nTokens < 5)  {
157 >                      
158 >    } else {
159 >
160 >        std::string atomTypeName = tokenizer.nextToken();
161 >        bool isDirectional = tokenizer.nextTokenAsBool();
162 >        bool isLJ = tokenizer.nextTokenAsBool();
163 >        bool isCharge = tokenizer.nextTokenAsBool();        
164 >        double mass = tokenizer.nextTokenAsDouble();
165 >        double epsilon;
166 >        double sigma;
167 >        double charge;
168 >        nTokens -= 5;
169 >
170 >        //parse epsilon and sigma
171 >        if (isLJ) {
172 >            if (nTokens >= 2) {
173 >                epsilon = tokenizer.nextTokenAsDouble();
174 >                sigma = tokenizer.nextTokenAsDouble();
175 >                nTokens -= 2;
176 >            } else {
177 >
178 >            }
179 >        }
180 >
181 >        //parse charge
182 >        if (isCharge) {
183 >            if (nTokens >= 1) {
184 >                charge = tokenizer.nextTokenAsDouble();
185 >                nTokens -= 1;
186 >            } else {
187 >
188 >            }
189 >        }
190 >        
191 >        AtomType* atomType;
192 >        if (isDirectional) {
193 >            atomType = new DirectionalAtomType();        
194 >        } else {
195 >            atomType = new AtomType();
196 >        }
197 >        
198 >        atomType->setName(atomTypeName);
199 >        atomType->setMass(mass);
200 >        
201 >        if (isLJ) {
202 >            atomType->setLennardJones();
203 >        }
204 >
205 >        if (isCharge) {
206 >            atomType->setCharge();
207 >        }
208 >        
209 >        atomType->setIdent(ident);
210 >        
211 >        atomType->complete();
212 >
213 >        int setLJStatus;
214 >        
215 >        //notify a new LJtype atom type is created
216 >        if (isLJ) {
217 >            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
218 >        }
219 >
220 >        int setChargeStatus;
221 >        if (isCharge) {
222 >            newChargeType(&ident, &charge, &setChargeStatus)
223 >        }
224 >
225 >        if (setLJStatus && setChargeStatus) {
226 >            //add atom type to AtomTypeContainer
227 >            addAtomType(atomTypeName, atomType);
228 >            ++ident;
229 >        } else {
230 >            //error in notifying fortran
231 >            delete atomType;
232 >        }
233 >    }    
234 >    
235   }
236  
237 < void DUFF::parseBondType(){
237 >
238 > void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
239 >    StringTokenizer tokenizer(line);
240 >    int nTokens = tokenizer.countTokens();  
241 >    
242 >    //in irectionalAtomTypeSection, a line at least contains 6 tokens
243 >    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
244 >    if (nTokens < 6) {
245 >
246 >    } else {
247 >
248 >
249 >        std::string atomTypeName = tokenizer.nextToken();
250 >        bool isDipole = tokenizer.nextTokenAsBool();
251 >        bool isSticky = tokenizer.nextTokenAsBool();
252 >        double Ixx = tokenizer.nextTokenAsDouble();
253 >        double Iyy = tokenizer.nextTokenAsDouble();
254 >        double Izz = tokenizer.nextTokenAsDouble();
255 >        nTokens -= 6;
256 >
257 >        AtomType* atomType = getAtomType(atomTypeName);
258 >        if (atomType == NULL) {
259 >
260 >        }
261 >        
262 >        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
263 >        if (dAtomType == NULL) {
264 >
265 >
266 >        }
267 >
268 >        if (isDipole) {
269 >            dAtomType->setDipole();
270 >        }
271 >
272 >        if (isSticky) {
273 >            dAtomType->setSticky();
274 >        }
275 >        
276 >        Mat3x3d inertialMat;
277 >        inertialMat(0, 0) = Ixx;
278 >        inertialMat(1, 1) = Ixx;
279 >        inertialMat(2, 2) = Ixx;        
280 >        dAtomType->setI(inertialMat);
281 >
282 >        //read dipole moment
283 >        double dipole;
284 >        if (isDipole) {
285 >            if (nTokens >= 1) {
286 >                dipole = tokenizer.nextTokenAsDouble();
287 >                nTokens -= 1;
288 >            } else {
289 >
290 >            }
291 >        }
292 >
293 >        //read sticky parameters
294 >        double w0;
295 >        double v0;
296 >        double v0p;
297 >        double rl;
298 >        double ru;
299 >        double rlp;
300 >        double rup;        
301 >        if (isSticky) {
302 >            if (nTokens >= 7) {
303 >                w0 = tokenizer.nextTokenAsDouble();
304 >                v0 = tokenizer.nextTokenAsDouble();
305 >                v0p = tokenizer.nextTokenAsDouble();
306 >                rl = tokenizer.nextTokenAsDouble();
307 >                ru = tokenizer.nextTokenAsDouble();
308 >                rlp = tokenizer.nextTokenAsDouble();
309 >                rup = tokenizer.nextTokenAsDouble();    
310 >                nTokens -= 7;
311 >            } else {
312 >
313 >            }
314 >        }
315 >
316 >
317 >        //notify fotran a newDipoleType is created
318 >        int ident = dAtomType->getIdent();
319 >        int setDipoleStatus;
320 >        if (isDipole) {
321 >            newDipoleType(&ident, &dipole, &setDipoleStatus);
322 >        }              
323 >
324 >        //notify fotran a StickyType is created
325 >        int setStickyStatus;
326 >        if (isSticky) {
327 >            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
328 >        }
329 >
330 >
331 >        if (!setDipoleStatus || !setStickyStatus) {
332 >
333 >        }
334 >        
335 >    }
336   }
337  
338 < void DUFF::parseBendType(){
338 > void DUFF::parseBondType(const std::string& line, int lineNo){
339 >
340 >    StringTokenizer tokenizer(line);
341 >    std::string at1;
342 >    std::string at2;
343 >    std::string bt;
344 >    BondType* bondType = NULL;
345 >    double b0;
346 >    
347 >    int nTokens = tokenizer.countTokens();
348 >
349 >    if (nTokens < 4) {
350 >
351 >        return;
352 >    }
353 >    
354 >    at1 = tokenizer.nextToken();
355 >    at2 = tokenizer.nextToken();
356 >    bt = tokenizer.nextToken();
357 >    b0 = tokenizer.nextTokenAsDouble();
358 >    nTokens -= 4;
359 >
360 >    //switch is a maintain nightmare
361 >    switch(bt) {
362 >        case "FixedBondType" :
363 >            bondType = new FixedBondType();
364 >            break;
365 >            
366 >        case "HarmonicBondType" :
367 >            if (nTokens < 1) {
368 >
369 >            } else {
370 >
371 >                double kb = tokenizer.nextTokenAsDouble();
372 >                bondType = new HarmonicBondType(b0, kb);
373 >            }
374 >
375 >            break;
376 >
377 >        case "CubicBondType" :
378 >            if (nTokens < 4) {
379 >
380 >            } else {
381 >
382 >                double k3 = tokenizer.nextTokenAsDouble();
383 >                double k2 = tokenizer.nextTokenAsDouble();
384 >                double k1 = tokenizer.nextTokenAsDouble();
385 >                double k0 = tokenizer.nextTokenAsDouble();
386 >                
387 >                bondType = new CubicBondType(b0, k3, k2, k1, k0);
388 >            }
389 >            break;
390 >            
391 >        case "QuadraticBondType" :
392 >            if (nTokens < 5) {
393 >
394 >            } else {
395 >
396 >                b0 = tokenizer.nextTokenAsDouble();
397 >                double k4 = tokenizer.nextTokenAsDouble();
398 >                double k3 = tokenizer.nextTokenAsDouble();
399 >                double k2 = tokenizer.nextTokenAsDouble();
400 >                double k1 = tokenizer.nextTokenAsDouble();
401 >                double k0 = tokenizer.nextTokenAsDouble();
402 >                
403 >                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
404 >            }
405 >            break;
406 >
407 >        case "PolynomialBondType " :
408 >            if (nTokens < 2 || nTokens % 2 != 0) {
409 >
410 >            } else {
411 >                int nPairs = nTokens / 2;
412 >                int power;
413 >                double coefficient;
414 >                PolynomialBondType pbt = new PolynomialBondType();
415 >                
416 >                for (int i = 0; i < nPairs; ++i) {
417 >                    power = tokenizer.nextTokenAsInt();
418 >                    coefficient = tokenizer.nextTokenAsDouble();
419 >                    pbt->setCoefficient(power, coefficient);
420 >                }
421 >            }
422 >            
423 >            break;
424 >
425 >        default:
426 >            
427 >    }
428 >
429 >    if (bondType != NULL) {
430 >        addBondType(at1, at2, bondType);
431 >    }
432   }
433  
434 < void DUFF::parseTorsionType(){
434 > void DUFF::parseBendType(const std::string& line, int lineNo){
435 >    StringTokenizer tokenizer(line);
436 >    std::string at1;
437 >    std::string at2;
438 >    std::string at3;
439 >    std::string bt;
440 >    double theta0;
441 >    BendType* bendType = NULL;
442  
443 +    int nTokens = tokenizer.countTokens();
444 +
445 +    if (nTokens < 5) {
446 +
447 +        return;
448 +    }
449 +    
450 +    at1 = tokenizer.nextToken();
451 +    at2 = tokenizer.nextToken();
452 +    at3 = tokenizer.nextToken();
453 +    bt = tokenizer.nextToken();
454 +    theta0 = tokenizer.nextTokenAsDouble();
455 +    nTokens -= 5;
456 +
457 +    //switch is a maintain nightmare
458 +    switch(bt) {
459 +            
460 +        case "HarmonicBendType" :
461 +            
462 +            if (nTokens < 1) {
463 +
464 +            } else {
465 +
466 +                double ktheta = tokenizer.nextTokenAsDouble();
467 +                bendType = new HarmonicBendType(theta0, ktheta);
468 +            }
469 +            break;
470 +        case "GhostBendType" :
471 +            if (nTokens < 1) {
472 +
473 +            } else {
474 +                double ktheta = tokenizer.nextTokenAsDouble();
475 +                bendType = new HarmonicBendType(theta0, ktheta);                
476 +            }
477 +            break;            
478 +
479 +        case "UreyBradleyBendType" :
480 +            if (nTokens < 3) {
481 +
482 +            } else {
483 +                double ktheta = tokenizer.nextTokenAsDouble();
484 +                double s0 =  tokenizer.nextTokenAsDouble();
485 +                double kub = tokenizer.nextTokenAsDouble();
486 +                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
487 +            }
488 +            break;
489 +            
490 +        case "CubicBendType" :
491 +            if (nTokens < 4) {
492 +
493 +            } else {
494 +
495 +                double k3 = tokenizer.nextTokenAsDouble();
496 +                double k2 = tokenizer.nextTokenAsDouble();
497 +                double k1 = tokenizer.nextTokenAsDouble();
498 +                double k0 = tokenizer.nextTokenAsDouble();
499 +                
500 +                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
501 +            }
502 +            break;
503 +            
504 +        case "QuadraticBendType" :
505 +            if (nTokens < 5) {
506 +
507 +            } else {
508 +
509 +                theta0 = tokenizer.nextTokenAsDouble();
510 +                double k4 = tokenizer.nextTokenAsDouble();
511 +                double k3 = tokenizer.nextTokenAsDouble();
512 +                double k2 = tokenizer.nextTokenAsDouble();
513 +                double k1 = tokenizer.nextTokenAsDouble();
514 +                double k0 = tokenizer.nextTokenAsDouble();
515 +                
516 +                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
517 +            }
518 +            break;
519 +
520 +        case "PolynomialBendType " :
521 +            if (nTokens < 2 || nTokens % 2 != 0) {
522 +
523 +            } else {
524 +                int nPairs = nTokens / 2;
525 +                int power;
526 +                double coefficient;
527 +                PolynomialBendType* pbt = new PolynomialBendType();
528 +                
529 +                for (int i = 0; i < nPairs; ++i) {
530 +                    power = tokenizer.nextTokenAsInt();
531 +                    coefficient = tokenizer.nextTokenAsDouble();
532 +                    pbt->setCoefficient(power, coefficient);
533 +                }
534 +            }
535 +            
536 +            break;
537 +
538 +        default:
539 +            
540 +    }
541 +
542 +    if (bendType != NULL) {
543 +        addBendType(at1, at2, at3, bendType);
544 +    }
545 +    
546   }
547  
548 < } //end namespace oopse
548 > void DUFF::parseTorsionType(const std::string& line, int lineNo){
549 >    StringTokenizer tokenizer(line);
550 >    std::string at1;
551 >    std::string at2;
552 >    std::string at3;
553 >    std::string at4;
554 >    std::string tt;
555 >    TorsionType* bendType = NULL;
556  
557 +    int nTokens = tokenizer.countTokens();
558 +
559 +    if (nTokens < 5) {
560 +
561 +        return;
562 +    }
563 +    
564 +    at1 = tokenizer.nextToken();
565 +    at2 = tokenizer.nextToken();
566 +    at3 = tokenizer.nextToken();
567 +    at4 = tokenizer.nextToken();
568 +    tt = tokenizer.nextToken();
569 +
570 +    nTokens -= 5;
571 +
572 +    switch(tt) {
573 +            
574 +        case "CubicTorsionType" :
575 +            if (nTokens < 4) {
576 +
577 +            } else {
578 +
579 +                double k3 = tokenizer.nextTokenAsDouble();
580 +                double k2 = tokenizer.nextTokenAsDouble();
581 +                double k1 = tokenizer.nextTokenAsDouble();
582 +                double k0 = tokenizer.nextTokenAsDouble();
583 +                
584 +                bendType = new CubicTorsionType(k3, k2, k1, k0);
585 +            }
586 +            break;
587 +            
588 +        case "QuadraticTorsionType" :
589 +            if (nTokens < 5) {
590 +
591 +            } else {
592 +
593 +                theta0 = tokenizer.nextTokenAsDouble();
594 +                double k4 = tokenizer.nextTokenAsDouble();
595 +                double k3 = tokenizer.nextTokenAsDouble();
596 +                double k2 = tokenizer.nextTokenAsDouble();
597 +                double k1 = tokenizer.nextTokenAsDouble();
598 +                double k0 = tokenizer.nextTokenAsDouble();
599 +                
600 +                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
601 +            }
602 +            break;
603 +
604 +        case "PolynomialTorsionType " :
605 +            if (nTokens < 2 || nTokens % 2 != 0) {
606 +
607 +            } else {
608 +                int nPairs = nTokens / 2;
609 +                int power;
610 +                double coefficient;
611 +                PolynomialTorsionType* pbt = new PolynomialTorsionType();
612 +                
613 +                for (int i = 0; i < nPairs; ++i) {
614 +                    power = tokenizer.nextTokenAsInt();
615 +                    coefficient = tokenizer.nextTokenAsDouble();
616 +                    pbt->setCoefficient(power, coefficient);
617 +                }
618 +            }
619 +            
620 +            break;
621 +        case "CharmmTorsionType" :
622 +            
623 +            if (nTokens < 3 || nTokens % 3 != 0) {
624 +
625 +            } else {
626 +                int nSets = nTokens / 3;
627 +  
628 +                CharmmTorsionType* ctt = new CharmmTorsionType();
629 +                
630 +                for (int i = 0; i < nSets; ++i) {
631 +                    double kchi = tokenizer.nextTokenAsDouble();
632 +                    int n = tokenizer.nextTokenAsInt();
633 +                    double delta = tokenizer.nextTokenAsDouble();
634 +    
635 +                    ctt->setCharmmTorsionParameter(kchi, n, delta);
636 +                }
637 +            }
638 +        default:
639 +            
640 +    }
641 +
642 +    if (bendType != NULL) {
643 +        addTorsionType(at1, at2, at3, bendType);
644 +    }
645 + }
646 +
647 + } //end namespace oopse

Diff Legend

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