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 1758 by tim, Fri Nov 19 17:56:32 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 >    ParseState currentSection = DUFF::UnknownSection;
80 >    
81 >    while(ffStream.getline(buffer, bufferSize)){
82          ++lineNo;
83  
84 <        //a line begins with '#' or '!' is comment which is skipped
85 <        if (line[0] != '#' ||line[0] != '!') {
86 <            StringTokenizer tokenizer(line);
87 <            
88 <            if (tokenizer.countTokens() >= 4) {
61 <                atomTypeName = tokenizer.nextToken();
62 <                mass = tokenizer.nextTokenAsDouble();
63 <                epsilon = tokenizer.nextTokenAsDouble();
64 <                sigma = tokenizer.nextTokenAsDouble();
84 >        line = trimSpaces(buffer);
85 >        //a line begins with "//" is comment
86 >        if ( line.empty() || (line.size() >= 2 && line[0] == '/' && line[1] == '/')) {
87 >            continue;
88 >        } else {
89  
90 <                atomType = new AtomType();
91 <                atomType->setName(atomTypeName);
92 <                atomType->setMass(mass);
93 <                
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();
90 >            switch(currentSection) {
91 >                case DUFF::AtomTypeSection :
92 >                    parseAtomType(line, lineNo);
93 >                    break;
94  
95 <                atomType->setIdent(ident);
96 <                atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
97 <                atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
77 <                atomType->complete();
78 <                //notify a new LJtype atom type is created
79 <                newLJtype(&ident, &sigma, &epsilon, &status);
95 >                case DUFF::BondTypeSection :
96 >                    parseBondType(line, lineNo);
97 >                    break;
98  
99 <                //add atom type to AtomTypeContainer
100 <                addAtomType(atomTypeName, atomType);
101 <                ++ident;
99 >                case DUFF::BendTypeSection :
100 >                    parseBendType(line, lineNo);
101 >                    break;
102 >
103 >                case DUFF::TorsionTypeSection :
104 >                    parseTorsionType(line, lineNo);
105 >                    break;
106 >                    
107 >                case DUFF::UnknownSection:
108 >                    StringTokenizer tokenizer(line);
109 >                    
110 >                    std::string keyword = tokenizer.nextToken();
111 >                    std::string section = tokenizer.nextToken();
112 >
113 >                    ParseState newSection = getSection(section);
114 >                    if (keyword != "begin" || keyword != "end") {
115 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
116 >                    } else if (keyword == "begin") {
117 >                        if (newSection == DUFF::UnknownSection) {
118 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
119 >                        } else {
120 >                            //enter a new section
121 >                            currentSection = newSection;
122 >                        }
123 >                        
124 >                    } else if (keyword == "end"){
125 >                        if (currentSection == newSection) ) {
126 >                            //leave a section
127 >                            currentSection = DUFF::UnknownSection;
128 >                        } else {                            
129 >                            std::cerr << "DUFF Parsing Error at line " << lineNo << ": " << line << std::endl;
130 >                        }
131 >
132 >                    }
133 >                    break;
134 >                default :
135                  
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();                
136              }
137              
138          }
# Line 99 | Line 142 | void DUFF::parseAtomType(){
142   }
143  
144  
145 < void DUFF::parseAtomType(){
145 > void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
146 >    StringTokenizer tokenizer(line);
147 >    int nTokens = tokenizer.countTokens();    
148 >
149 >    //in AtomTypeSection, a line at least contains 5 tokens
150 >    //atomTypeName, is Directional, isLJ, isCharge and mass
151 >    if (nTokens < 5)  {
152 >        sprintf( painCave.errMsg,
153 >               "Not enough tokens when parsing DUFF Force Field : %s\n"
154 >               "in line %d : %s\n",
155 >               filename.c_str(), lineNo, line);
156 >        painCave.severity = OOPSE_ERROR;
157 >        painCave.isFatal = 1;
158 >        simError();                
159 >
160 >        
161 >    } else {
162 >
163 >        std::string atomTypeName = tokenizer.nextToken();
164 >        bool isDirectional = tokenizer.nextTokenAsBool();
165 >        bool isLJ = tokenizer.nextTokenAsBool();
166 >        bool isCharge = tokenizer.nextTokenAsBool();        
167 >        double mass = tokenizer.nextTokenAsDouble();
168 >        double epsilon;
169 >        double sigma;
170 >        double charge;
171 >        nTokens -= 5;
172 >
173 >        //parse epsilon and sigma
174 >        if (isLJ) {
175 >            if (nTokens >= 2) {
176 >                epsilon = tokenizer.nextTokenAsDouble();
177 >                sigma = tokenizer.nextTokenAsDouble();
178 >                nTokens -= 2;
179 >            } else {
180 >
181 >            }
182 >        }
183 >
184 >        //parse charge
185 >        if (isCharge) {
186 >            if (nTokens >= 1) {
187 >                charge = tokenizer.nextTokenAsDouble();
188 >                nTokens -= 1;
189 >            } else {
190 >
191 >            }
192 >        }
193 >        
194 >        AtomType* atomType;
195 >        if (isDirectional) {
196 >            atomType = new DirectionalAtomType();        
197 >        } else {
198 >            atomType = new AtomType();
199 >        }
200 >        
201 >        atomType->setName(atomTypeName);
202 >        atomType->setMass(mass);
203 >        
204 >        //by default, all of the properties in AtomTypeProperties is set to 0
205 >        //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
206 >        atomType->setLennardJones();
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 irectionalAtomTypeSection, a line at least contains 6 tokens
242 >    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
243 >    if (nTokens < 6) {
244 >
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 >        dAtomType->setDipole(isDipole);
268 >        dAtomType->setSticky(isSticky);
269 >
270 >        Mat3x3d inertialMat;
271 >        inertialMat(0, 0) = Ixx;
272 >        inertialMat(1, 1) = Ixx;
273 >        inertialMat(2, 2) = Ixx;        
274 >        dAtomType->setI(inertialMat);
275 >
276 >        //read dipole moment
277 >        double dipole;
278 >        if (isDipole) {
279 >            if (nTokens >= 1) {
280 >                dipole = tokenizer.nextTokenAsDouble();
281 >                nTokens -= 1;
282 >            } else {
283 >
284 >            }
285 >        }
286 >
287 >        //read sticky parameters
288 >        double w0;
289 >        double v0;
290 >        double v0p;
291 >        double rl;
292 >        double ru;
293 >        double rlp;
294 >        double rup;        
295 >        if (isSticky) {
296 >            if (nTokens >= 7) {
297 >                w0 = tokenizer.nextTokenAsDouble();
298 >                v0 = tokenizer.nextTokenAsDouble();
299 >                v0p = tokenizer.nextTokenAsDouble();
300 >                rl = tokenizer.nextTokenAsDouble();
301 >                ru = tokenizer.nextTokenAsDouble();
302 >                rlp = tokenizer.nextTokenAsDouble();
303 >                rup = tokenizer.nextTokenAsDouble();    
304 >                nTokens -= 7;
305 >            } else {
306 >
307 >            }
308 >        }
309 >
310 >
311 >        //notify fotran a newDipoleType is created
312 >        int ident = dAtomType->getIdent();
313 >        int setDipoleStatus;
314 >        if (isDipole) {
315 >            newDipoleType(&ident, &dipole, &setDipoleStatus);
316 >        }              
317 >
318 >        //notify fotran a StickyType is created
319 >        int setStickyStatus;
320 >        if (isSticky) {
321 >            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
322 >        }
323 >
324 >
325 >        if (!setDipoleStatus || !setStickyStatus) {
326 >
327 >        }
328 >        
329 >    }
330   }
331  
332 < void DUFF::parseBendType(){
332 > void DUFF::parseBondType(const std::string& line, int lineNo){
333 >
334 >    StringTokenizer tokenizer(line);
335 >    std::string at1;
336 >    std::string at2;
337 >    std::string bt;
338 >    BondType* bondType = NULL;
339 >    double b0;
340 >    
341 >    int nTokens = tokenizer.countTokens();
342 >
343 >    if (nTokens < 4) {
344 >
345 >        return;
346 >    }
347 >    
348 >    at1 = tokenizer.nextToken();
349 >    at2 = tokenizer.nextToken();
350 >    bt = tokenizer.nextToken();
351 >    b0 = tokenizer.nextTokenAsDouble();
352 >    nTokens -= 4;
353 >
354 >    //switch is a maintain nightmare
355 >    switch(bt) {
356 >        case "FixedBondType" :
357 >            bondType = new FixedBondType();
358 >            break;
359 >            
360 >        case "HarmonicBondType" :
361 >            if (nTokens < 1) {
362 >
363 >            } else {
364 >
365 >                double kb = tokenizer.nextTokenAsDouble();
366 >                bondType = new HarmonicBondType(b0, kb);
367 >            }
368 >
369 >            break;
370 >
371 >        case "CubicBondType" :
372 >            if (nTokens < 4) {
373 >
374 >            } else {
375 >
376 >                double k3 = tokenizer.nextTokenAsDouble();
377 >                double k2 = tokenizer.nextTokenAsDouble();
378 >                double k1 = tokenizer.nextTokenAsDouble();
379 >                double k0 = tokenizer.nextTokenAsDouble();
380 >                
381 >                bondType = new CubicBondType(b0, k3, k2, k1, k0);
382 >            }
383 >            break;
384 >            
385 >        case "QuadraticBondType" :
386 >            if (nTokens < 5) {
387 >
388 >            } else {
389 >
390 >                b0 = tokenizer.nextTokenAsDouble();
391 >                double k4 = tokenizer.nextTokenAsDouble();
392 >                double k3 = tokenizer.nextTokenAsDouble();
393 >                double k2 = tokenizer.nextTokenAsDouble();
394 >                double k1 = tokenizer.nextTokenAsDouble();
395 >                double k0 = tokenizer.nextTokenAsDouble();
396 >                
397 >                bondType = new QuadraticBondType(b0, k4, k3, k2, k1, k0);
398 >            }
399 >            break;
400 >
401 >        case "PolynomialBondType " :
402 >            if (nTokens < 2 || nTokens % 2 != 0) {
403 >
404 >            } else {
405 >                int nPairs = nTokens / 2;
406 >                int power;
407 >                double coefficient;
408 >                PolynomialBondType pbt = new PolynomialBondType();
409 >                
410 >                for (int i = 0; i < nPairs; ++i) {
411 >                    power = tokenizer.nextTokenAsInt();
412 >                    coefficient = tokenizer.nextTokenAsDouble();
413 >                    pbt->setCoefficient(power, coefficient);
414 >                }
415 >            }
416 >            
417 >            break;
418 >
419 >        default:
420 >            
421 >    }
422 >
423 >    if (bondType != NULL) {
424 >        addBondType(at1, at2, bondType);
425 >    }
426   }
427  
428 < void DUFF::parseTorsionType(){
428 > void DUFF::parseBendType(const std::string& line, int lineNo){
429 >    StringTokenizer tokenizer(line);
430 >    std::string at1;
431 >    std::string at2;
432 >    std::string at3;
433 >    std::string bt;
434 >    double theta0;
435 >    BendType* bendType = NULL;
436  
437 +    int nTokens = tokenizer.countTokens();
438 +
439 +    if (nTokens < 5) {
440 +
441 +        return;
442 +    }
443 +    
444 +    at1 = tokenizer.nextToken();
445 +    at2 = tokenizer.nextToken();
446 +    at3 = tokenizer.nextToken();
447 +    bt = tokenizer.nextToken();
448 +    theta0 = tokenizer.nextTokenAsDouble();
449 +    nTokens -= 5;
450 +
451 +    //switch is a maintain nightmare
452 +    switch(bt) {
453 +            
454 +        case "HarmonicBendType" :
455 +            
456 +            if (nTokens < 1) {
457 +
458 +            } else {
459 +
460 +                double ktheta = tokenizer.nextTokenAsDouble();
461 +                bendType = new HarmonicBendType(theta0, ktheta);
462 +            }
463 +            break;
464 +        case "GhostBendType" :
465 +            if (nTokens < 1) {
466 +
467 +            } else {
468 +                double ktheta = tokenizer.nextTokenAsDouble();
469 +                bendType = new HarmonicBendType(theta0, ktheta);                
470 +            }
471 +            break;            
472 +
473 +        case "UreyBradleyBendType" :
474 +            if (nTokens < 3) {
475 +
476 +            } else {
477 +                double ktheta = tokenizer.nextTokenAsDouble();
478 +                double s0 =  tokenizer.nextTokenAsDouble();
479 +                double kub = tokenizer.nextTokenAsDouble();
480 +                bendType = new UreyBradleyBendType(theta0, ktheta, s0, kub);                
481 +            }
482 +            break;
483 +            
484 +        case "CubicBendType" :
485 +            if (nTokens < 4) {
486 +
487 +            } else {
488 +
489 +                double k3 = tokenizer.nextTokenAsDouble();
490 +                double k2 = tokenizer.nextTokenAsDouble();
491 +                double k1 = tokenizer.nextTokenAsDouble();
492 +                double k0 = tokenizer.nextTokenAsDouble();
493 +                
494 +                bendType = new CubicBendType(theta0, k3, k2, k1, k0);
495 +            }
496 +            break;
497 +            
498 +        case "QuadraticBendType" :
499 +            if (nTokens < 5) {
500 +
501 +            } else {
502 +
503 +                theta0 = tokenizer.nextTokenAsDouble();
504 +                double k4 = tokenizer.nextTokenAsDouble();
505 +                double k3 = tokenizer.nextTokenAsDouble();
506 +                double k2 = tokenizer.nextTokenAsDouble();
507 +                double k1 = tokenizer.nextTokenAsDouble();
508 +                double k0 = tokenizer.nextTokenAsDouble();
509 +                
510 +                bendType = new QuadraticBendType(theta0, k4, k3, k2, k1, k0);
511 +            }
512 +            break;
513 +
514 +        case "PolynomialBendType " :
515 +            if (nTokens < 2 || nTokens % 2 != 0) {
516 +
517 +            } else {
518 +                int nPairs = nTokens / 2;
519 +                int power;
520 +                double coefficient;
521 +                PolynomialBendType* pbt = new PolynomialBendType();
522 +                
523 +                for (int i = 0; i < nPairs; ++i) {
524 +                    power = tokenizer.nextTokenAsInt();
525 +                    coefficient = tokenizer.nextTokenAsDouble();
526 +                    pbt->setCoefficient(power, coefficient);
527 +                }
528 +            }
529 +            
530 +            break;
531 +
532 +        default:
533 +            
534 +    }
535 +
536 +    if (bendType != NULL) {
537 +        addBendType(at1, at2, at3, bendType);
538 +    }
539 +    
540   }
541  
542 < } //end namespace oopse
542 > void DUFF::parseTorsionType(const std::string& line, int lineNo){
543 >    StringTokenizer tokenizer(line);
544 >    std::string at1;
545 >    std::string at2;
546 >    std::string at3;
547 >    std::string at4;
548 >    std::string tt;
549 >    TorsionType* bendType = NULL;
550  
551 +    int nTokens = tokenizer.countTokens();
552 +
553 +    if (nTokens < 5) {
554 +
555 +        return;
556 +    }
557 +    
558 +    at1 = tokenizer.nextToken();
559 +    at2 = tokenizer.nextToken();
560 +    at3 = tokenizer.nextToken();
561 +    at4 = tokenizer.nextToken();
562 +    tt = tokenizer.nextToken();
563 +
564 +    nTokens -= 5;
565 +
566 +    switch(tt) {
567 +            
568 +        case "CubicTorsionType" :
569 +            if (nTokens < 4) {
570 +
571 +            } else {
572 +
573 +                double k3 = tokenizer.nextTokenAsDouble();
574 +                double k2 = tokenizer.nextTokenAsDouble();
575 +                double k1 = tokenizer.nextTokenAsDouble();
576 +                double k0 = tokenizer.nextTokenAsDouble();
577 +                
578 +                bendType = new CubicTorsionType(k3, k2, k1, k0);
579 +            }
580 +            break;
581 +            
582 +        case "QuadraticTorsionType" :
583 +            if (nTokens < 5) {
584 +
585 +            } else {
586 +
587 +                theta0 = tokenizer.nextTokenAsDouble();
588 +                double k4 = tokenizer.nextTokenAsDouble();
589 +                double k3 = tokenizer.nextTokenAsDouble();
590 +                double k2 = tokenizer.nextTokenAsDouble();
591 +                double k1 = tokenizer.nextTokenAsDouble();
592 +                double k0 = tokenizer.nextTokenAsDouble();
593 +                
594 +                bendType = new QuadraticTorsionType( k4, k3, k2, k1, k0);
595 +            }
596 +            break;
597 +
598 +        case "PolynomialTorsionType " :
599 +            if (nTokens < 2 || nTokens % 2 != 0) {
600 +
601 +            } else {
602 +                int nPairs = nTokens / 2;
603 +                int power;
604 +                double coefficient;
605 +                PolynomialTorsionType* pbt = new PolynomialTorsionType();
606 +                
607 +                for (int i = 0; i < nPairs; ++i) {
608 +                    power = tokenizer.nextTokenAsInt();
609 +                    coefficient = tokenizer.nextTokenAsDouble();
610 +                    pbt->setCoefficient(power, coefficient);
611 +                }
612 +            }
613 +            
614 +            break;
615 +        case "CharmmTorsionType" :
616 +            
617 +            if (nTokens < 3 || nTokens % 3 != 0) {
618 +
619 +            } else {
620 +                int nSets = nTokens / 3;
621 +  
622 +                CharmmTorsionType* ctt = new CharmmTorsionType();
623 +                
624 +                for (int i = 0; i < nSets; ++i) {
625 +                    double kchi = tokenizer.nextTokenAsDouble();
626 +                    int n = tokenizer.nextTokenAsInt();
627 +                    double delta = tokenizer.nextTokenAsDouble();
628 +    
629 +                    ctt->setCharmmTorsionParameter(kchi, n, delta);
630 +                }
631 +            }
632 +        default:
633 +            
634 +    }
635 +
636 +    if (bendType != NULL) {
637 +        addTorsionType(at1, at2, at3, bendType);
638 +    }
639 + }
640 +
641 + } //end namespace oopse

Diff Legend

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