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

Diff Legend

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