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 1757 by tim, Fri Nov 19 00:11:33 2004 UTC vs.
Revision 1770 by tim, Tue Nov 23 17:53:43 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 + DUFF::DUFF(){
43  
44 +    //the order of adding section parsers are important
45 +    //DirectionalAtomTypesSectionParser should be added before AtomTypesSectionParser Since
46 +    //These two section parsers will actually create "real" AtomTypes (AtomTypesSectionParser will create
47 +    //AtomType and DirectionalAtomTypesSectionParser will creat DirectionalAtomType which is a subclass
48 +    //of AtomType, therefore it should come first). Other AtomTypes Section Parser will not create the
49 +    //"real" AtomType, they only add and set some attribute of the AtomType. Thus their order are not
50 +    //important. AtomTypesSectionParser should be added before other atom type section parsers.
51 +    //Make sure they are added after DirectionalAtomTypesSectionParser and AtomTypesSectionParser.
52 +    //The order of BondTypesSectionParser, BendTypesSectionParser and TorsionTypesSectionParser are
53 +    //not important.
54 +    spMan_.push_back(new DirectionalAtomTypesSectionParser());
55 +    spMan_.push_back(new AtomTypesSectionParser());
56 +    spMan_.push_back(new LennardJonesAtomTypesSectionParser());
57 +    spMan_.push_back(new ElectrostaticAtomTypesSectionParser());
58 +    spMan_.push_back(new EAMAtomTypesSectionParser());
59 +    spMan_.push_back(new StickyAtomTypesSectionParser());
60 +    spMan_.push_back(new BondTypesSectionParser());
61 +    spMan_.push_back(new BendTypesSectionParser());
62 +    spMan_.push_back(new TorsionTypesSectionParser());
63 +    
64 + }
65 +
66 + void DUFF::parse(const std::string& filename) {
67 +    ifstrstream* ffStream;
68 +    ffStream = openForceFiledFile(filename);
69 +
70 +    spMan_.parse(*ffStream, *this);
71 +
72 +    typename ForceField::AtomTypeContainer::MapTypeIterator i;
73 +    AtomType* at;
74 +
75 +    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
76 +        at->makeFortranAtomType();
77 +    }
78 +
79 +    for (at = atomTypeCont_.beginType(i); at != NULL; at = atomTypeCont_.nextType(i)) {
80 +        at->complete();
81 +    }
82 +    
83 + }
84 +
85 + /*
86   ParseState DUFF::getSection(const std::string& section) {
87      ParseState result;
88      
# Line 44 | Line 90 | ParseState DUFF::getSection(const std::string& section
90          case "AtomTypes" :
91              result = DUFF::AtomTypeSection;
92              break;
93 +        case "DirectionalAtomTypes" :
94 +            result = DUFF::DirectionalAtomTypeSection;
95 +            break;
96  
97          case "BondTypes" :
98              result = DUFF::BondTypeSection;
# Line 70 | Line 119 | void DUFF::parse(const std::string& filename) {
119      std::string line;
120      char buffer[bufferSize];
121      int lineNo = 0;
122 +    int atomIdent = getNAtomType() + 1;  //atom's indent begins from 1 (since only fortran's array begins from 1)
123      ParseState currentSection = DUFF::UnknownSection;
124      
125      while(ffStream.getline(buffer, bufferSize)){
# Line 83 | Line 133 | void DUFF::parse(const std::string& filename) {
133  
134              switch(currentSection) {
135                  case DUFF::AtomTypeSection :
136 <                    parseAtomType(line, lineNo);
136 >                    parseAtomType(line, lineNo, atomIdent);
137                      break;
138  
139 +                case DUFF::DirectionalAtomTypeSection :
140 +                    parseDirectionalAtomType(line, lineNo);
141 +                    break;
142 +                    
143                  case DUFF::BondTypeSection :
144                      parseBondType(line, lineNo);
145                      break;
# Line 135 | Line 189 | void DUFF::parse(const std::string& filename) {
189      delete ffStream;
190   }
191  
192 <
139 < void DUFF::parseAtomType(const std::string& line, int lineNo){
192 > void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){
193      StringTokenizer tokenizer(line);
194 <    
142 <    if (tokenizer.countTokens() >= 4) {
143 <        atomTypeName = tokenizer.nextToken();
144 <        mass = tokenizer.nextTokenAsFloat();
145 <        epsilon = tokenizer.nextTokenAsFloat();
146 <        sigma = tokenizer.nextTokenAsFloat();
194 >    int nTokens = tokenizer.countTokens();    
195  
196 <        atomType = new AtomType();
196 >    //in AtomTypeSection, a line at least contains 5 tokens
197 >    //atomTypeName, is Directional, isLJ, isCharge and mass
198 >    if (nTokens < 5)  {
199 >                      
200 >    } else {
201 >
202 >        std::string atomTypeName = tokenizer.nextToken();
203 >        bool isDirectional = tokenizer.nextTokenAsBool();
204 >        bool isLJ = tokenizer.nextTokenAsBool();
205 >        bool isCharge = tokenizer.nextTokenAsBool();        
206 >        double mass = tokenizer.nextTokenAsDouble();
207 >        double epsilon;
208 >        double sigma;
209 >        double charge;
210 >        nTokens -= 5;
211 >
212 >        //parse epsilon and sigma
213 >        if (isLJ) {
214 >            if (nTokens >= 2) {
215 >                epsilon = tokenizer.nextTokenAsDouble();
216 >                sigma = tokenizer.nextTokenAsDouble();
217 >                nTokens -= 2;
218 >            } else {
219 >
220 >            }
221 >        }
222 >
223 >        //parse charge
224 >        if (isCharge) {
225 >            if (nTokens >= 1) {
226 >                charge = tokenizer.nextTokenAsDouble();
227 >                nTokens -= 1;
228 >            } else {
229 >
230 >            }
231 >        }
232 >        
233 >        AtomType* atomType;
234 >        if (isDirectional) {
235 >            atomType = new DirectionalAtomType();        
236 >        } else {
237 >            atomType = new AtomType();
238 >        }
239 >        
240          atomType->setName(atomTypeName);
241          atomType->setMass(mass);
242          
243 <        //by default, all of the properties in AtomTypeProperties is set to 0
244 <        //In  Lennard-Jones Force Field, we only need to set Lennard-Jones to true
245 <        atomType->setLennardJones();
243 >        if (isLJ) {
244 >            atomType->setLennardJones();
245 >        }
246  
247 <        atomType->setIdent(ident);
248 <        atomType->addProperty(new DoubleGenericData("Epsilon", epsilon));                                
249 <        atomType->addProperty(new DoubleGenericData("Sigma", sigma));                                
247 >        if (isCharge) {
248 >            atomType->setCharge();
249 >        }
250 >        
251 >        atomType->setIdent(ident);
252 >        
253          atomType->complete();
160 –        //notify a new LJtype atom type is created
161 –        newLJtype(&ident, &sigma, &epsilon, &status);
254  
255 <        //add atom type to AtomTypeContainer
164 <        addAtomType(atomTypeName, atomType);
165 <        ++ident;
255 >        int setLJStatus;
256          
257 <    } else {
258 <        sprintf( painCave.errMsg,
259 <               "Not enough tokens when parsing Lennard-Jones Force Field : %s\n"
260 <               "in line %d : %s\n",
261 <               filename.c_str(), lineNo, line);
262 <        painCave.severity = OOPSE_ERROR;
263 <        painCave.isFatal = 1;
264 <        simError();                
257 >        //notify a new LJtype atom type is created
258 >        if (isLJ) {
259 >            newLJtype(&ident, &sigma, &epsilon, &setLJStatus);
260 >        }
261 >
262 >        int setChargeStatus;
263 >        if (isCharge) {
264 >            newChargeType(&ident, &charge, &setChargeStatus)
265 >        }
266 >
267 >        if (setLJStatus && setChargeStatus) {
268 >            //add atom type to AtomTypeContainer
269 >            addAtomType(atomTypeName, atomType);
270 >            ++ident;
271 >        } else {
272 >            //error in notifying fortran
273 >            delete atomType;
274 >        }
275      }    
276 +    
277   }
278  
279 +
280 + void DUFF::parseDirectionalAtomType(const std::string& line, int lineNo) {
281 +    StringTokenizer tokenizer(line);
282 +    int nTokens = tokenizer.countTokens();  
283 +    
284 +    //in DirectionalAtomTypeSection, a line at least contains 6 tokens
285 +    //AtomTypeName, isDipole, isSticky, I_xx, I_yy and I_zz
286 +    if (nTokens < 6) {
287 +        std::cerr << "Not enought tokens" << std::endl;
288 +    } else {
289 +
290 +
291 +        std::string atomTypeName = tokenizer.nextToken();
292 +        bool isDipole = tokenizer.nextTokenAsBool();
293 +        bool isSticky = tokenizer.nextTokenAsBool();
294 +        double Ixx = tokenizer.nextTokenAsDouble();
295 +        double Iyy = tokenizer.nextTokenAsDouble();
296 +        double Izz = tokenizer.nextTokenAsDouble();
297 +        nTokens -= 6;
298 +
299 +        AtomType* atomType = getAtomType(atomTypeName);
300 +        if (atomType == NULL) {
301 +
302 +        }
303 +        
304 +        DirectionalAtomType* dAtomType = dynamic_cast<DirectionalAtomType*>(atomType);
305 +        if (dAtomType == NULL) {
306 +
307 +
308 +        }
309 +
310 +        if (isDipole) {
311 +            dAtomType->setDipole();
312 +        }
313 +
314 +        if (isSticky) {
315 +            dAtomType->setSticky();
316 +        }
317 +        
318 +        Mat3x3d inertialMat;
319 +        inertialMat(0, 0) = Ixx;
320 +        inertialMat(1, 1) = Iyy;
321 +        inertialMat(2, 2) = Izz;        
322 +        dAtomType->setI(inertialMat);
323 +
324 +        //read dipole moment
325 +        double dipole;
326 +        if (isDipole) {
327 +            if (nTokens >= 1) {
328 +                dipole = tokenizer.nextTokenAsDouble();
329 +                nTokens -= 1;
330 +            } else {
331 +
332 +            }
333 +        }
334 +
335 +        //read sticky parameters
336 +        double w0;
337 +        double v0;
338 +        double v0p;
339 +        double rl;
340 +        double ru;
341 +        double rlp;
342 +        double rup;        
343 +        if (isSticky) {
344 +            if (nTokens >= 7) {
345 +                w0 = tokenizer.nextTokenAsDouble();
346 +                v0 = tokenizer.nextTokenAsDouble();
347 +                v0p = tokenizer.nextTokenAsDouble();
348 +                rl = tokenizer.nextTokenAsDouble();
349 +                ru = tokenizer.nextTokenAsDouble();
350 +                rlp = tokenizer.nextTokenAsDouble();
351 +                rup = tokenizer.nextTokenAsDouble();    
352 +                nTokens -= 7;
353 +            } else {
354 +
355 +            }
356 +        }
357 +
358 +
359 +        //notify fotran a newDipoleType is created
360 +        int ident = dAtomType->getIdent();
361 +        int setDipoleStatus;
362 +        if (isDipole) {
363 +            newDipoleType(&ident, &dipole, &setDipoleStatus);
364 +        }              
365 +
366 +        //notify fotran a StickyType is created
367 +        int setStickyStatus;
368 +        if (isSticky) {
369 +            makeStickyType( &w0, &v0, &v0p, &rl, &ru, &rlp, &rup);
370 +        }
371 +
372 +
373 +        if (!setDipoleStatus || !setStickyStatus) {
374 +
375 +        }
376 +        
377 +    }
378 + }
379 +
380   void DUFF::parseBondType(const std::string& line, int lineNo){
381  
382      StringTokenizer tokenizer(line);
# Line 199 | Line 401 | void DUFF::parseBondType(const std::string& line, int
401  
402      //switch is a maintain nightmare
403      switch(bt) {
404 <        case "FixedBondType" :
404 >        case "Fixed" :
405              bondType = new FixedBondType();
406              break;
407              
408 <        case "HarmonicBondType" :
408 >        case "Harmonic" :
409              if (nTokens < 1) {
410  
411              } else {
# Line 214 | Line 416 | void DUFF::parseBondType(const std::string& line, int
416  
417              break;
418  
419 <        case "CubicBondType" :
419 >        case "Cubic" :
420              if (nTokens < 4) {
421  
422              } else {
# Line 228 | Line 430 | void DUFF::parseBondType(const std::string& line, int
430              }
431              break;
432              
433 <        case "QuadraticBondType" :
433 >        case "Quartic" :
434              if (nTokens < 5) {
435  
436              } else {
# Line 244 | Line 446 | void DUFF::parseBondType(const std::string& line, int
446              }
447              break;
448  
449 <        case "PolynomialBondType " :
449 >        case "Polynomial" :
450              if (nTokens < 2 || nTokens % 2 != 0) {
451  
452              } else {
# Line 256 | Line 458 | void DUFF::parseBondType(const std::string& line, int
458                  for (int i = 0; i < nPairs; ++i) {
459                      power = tokenizer.nextTokenAsInt();
460                      coefficient = tokenizer.nextTokenAsDouble();
461 <
260 <                    if (power < 0) {
261 <
262 <                    } else {
263 <                        pbt->setCoefficient(power, coefficient);
264 <                    }
461 >                    pbt->setCoefficient(power, coefficient);
462                  }
463              }
464              
# Line 302 | Line 499 | void DUFF::parseBendType(const std::string& line, int
499      //switch is a maintain nightmare
500      switch(bt) {
501              
502 <        case "HarmonicBendType" :
502 >        case "Harmonic" :
503              
504              if (nTokens < 1) {
505  
# Line 312 | Line 509 | void DUFF::parseBendType(const std::string& line, int
509                  bendType = new HarmonicBendType(theta0, ktheta);
510              }
511              break;
512 <        case "GhostBendType" :
512 >        case "GhostBend" :
513              if (nTokens < 1) {
514  
515              } else {
# Line 321 | Line 518 | void DUFF::parseBendType(const std::string& line, int
518              }
519              break;            
520  
521 <        case "UreyBradleyBendType" :
521 >        case "UreyBradley" :
522              if (nTokens < 3) {
523  
524              } else {
# Line 332 | Line 529 | void DUFF::parseBendType(const std::string& line, int
529              }
530              break;
531              
532 <        case "CubicBendType" :
532 >        case "Cubic" :
533              if (nTokens < 4) {
534  
535              } else {
# Line 346 | Line 543 | void DUFF::parseBendType(const std::string& line, int
543              }
544              break;
545              
546 <        case "QuadraticBendType" :
546 >        case "Quartic" :
547              if (nTokens < 5) {
548  
549              } else {
# Line 362 | Line 559 | void DUFF::parseBendType(const std::string& line, int
559              }
560              break;
561  
562 <        case "PolynomialBendType " :
562 >        case "Polynomial" :
563              if (nTokens < 2 || nTokens % 2 != 0) {
564  
565              } else {
# Line 374 | Line 571 | void DUFF::parseBendType(const std::string& line, int
571                  for (int i = 0; i < nPairs; ++i) {
572                      power = tokenizer.nextTokenAsInt();
573                      coefficient = tokenizer.nextTokenAsDouble();
574 <
378 <                    if (power < 0) {
379 <
380 <                    } else {
381 <                        pbt->setCoefficient(power, coefficient);
382 <                    }
574 >                    pbt->setCoefficient(power, coefficient);
575                  }
576              }
577              
# Line 402 | Line 594 | void DUFF::parseTorsionType(const std::string& line, i
594      std::string at3;
595      std::string at4;
596      std::string tt;
597 <    TorsionType* bendType = NULL;
597 >    TorsionType* torsionType = NULL;
598  
599      int nTokens = tokenizer.countTokens();
600  
# Line 421 | Line 613 | void DUFF::parseTorsionType(const std::string& line, i
613  
614      switch(tt) {
615              
616 <        case "CubicTorsionType" :
616 >        case "Cubic" :
617              if (nTokens < 4) {
618  
619              } else {
# Line 435 | Line 627 | void DUFF::parseTorsionType(const std::string& line, i
627              }
628              break;
629              
630 <        case "QuadraticTorsionType" :
630 >        case "Quartic" :
631              if (nTokens < 5) {
632  
633              } else {
# Line 451 | Line 643 | void DUFF::parseTorsionType(const std::string& line, i
643              }
644              break;
645  
646 <        case "PolynomialTorsionType " :
646 >        case "Polynomial" :
647              if (nTokens < 2 || nTokens % 2 != 0) {
648  
649              } else {
# Line 463 | Line 655 | void DUFF::parseTorsionType(const std::string& line, i
655                  for (int i = 0; i < nPairs; ++i) {
656                      power = tokenizer.nextTokenAsInt();
657                      coefficient = tokenizer.nextTokenAsDouble();
658 <
467 <                    if (power < 0) {
468 <
469 <                    } else {
470 <                        pbt->setCoefficient(power, coefficient);
471 <                    }
658 >                    pbt->setCoefficient(power, coefficient);
659                  }
660              }
661              
662              break;
663 <        case "CharmmTorsionType" :
663 >        case "Charmm" :
664              
665              if (nTokens < 3 || nTokens % 3 != 0) {
666  
667              } else {
668                  int nSets = nTokens / 3;
669    
483 –                double kchi;
484 –                int n;
485 –                double delta;                
670                  CharmmTorsionType* ctt = new CharmmTorsionType();
671                  
672                  for (int i = 0; i < nSets; ++i) {
673 <                    kchi = tokenizer.nextTokenAsDouble();
674 <                    n = tokenizer.nextTokenAsInt();
675 <                    delta = tokenizer.nextTokenAsDouble();
673 >                    double kchi = tokenizer.nextTokenAsDouble();
674 >                    int n = tokenizer.nextTokenAsInt();
675 >                    double delta = tokenizer.nextTokenAsDouble();
676      
677 <                    if (n < 0 || n > 6) {
494 <
495 <                    } else {
496 <                        ctt->setCharmmTorsionParameter(kchi, n, delta);
497 <                    }
677 >                    ctt->setCharmmTorsionParameter(kchi, n, delta);
678                  }
679              }
680          default:
681              
682      }
683  
684 <    if (bendType != NULL) {
685 <        addTorsionType(at1, at2, at3, bendType);
684 >    if (torsionType != NULL) {
685 >        addTorsionType(at1, at2, at3, at4, torsionType);
686      }
687   }
688 <
688 > */
689   } //end namespace oopse
510 –

Diff Legend

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