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 |
|
|
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)){ |
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; |
189 |
|
delete ffStream; |
190 |
|
} |
191 |
|
|
144 |
– |
|
192 |
|
void DUFF::parseAtomType(const std::string& line, int lineNo, int& ident){ |
193 |
|
StringTokenizer tokenizer(line); |
194 |
|
int nTokens = tokenizer.countTokens(); |
196 |
|
//in AtomTypeSection, a line at least contains 5 tokens |
197 |
|
//atomTypeName, is Directional, isLJ, isCharge and mass |
198 |
|
if (nTokens < 5) { |
199 |
< |
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 |
< |
|
199 |
> |
|
200 |
|
} else { |
201 |
|
|
202 |
|
std::string atomTypeName = tokenizer.nextToken(); |
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 |
+ |
if (isCharge) { |
248 |
+ |
atomType->setCharge(); |
249 |
+ |
} |
250 |
+ |
|
251 |
|
atomType->setIdent(ident); |
252 |
|
|
253 |
|
atomType->complete(); |
281 |
|
StringTokenizer tokenizer(line); |
282 |
|
int nTokens = tokenizer.countTokens(); |
283 |
|
|
284 |
< |
//in irectionalAtomTypeSection, a line at least contains 6 tokens |
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 |
< |
|
287 |
> |
std::cerr << "Not enought tokens" << std::endl; |
288 |
|
} else { |
289 |
|
|
290 |
|
|
307 |
|
|
308 |
|
} |
309 |
|
|
310 |
< |
dAtomType->setDipole(isDipole); |
311 |
< |
dAtomType->setSticky(isSticky); |
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) = Ixx; |
321 |
< |
inertialMat(2, 2) = Ixx; |
320 |
> |
inertialMat(1, 1) = Iyy; |
321 |
> |
inertialMat(2, 2) = Izz; |
322 |
|
dAtomType->setI(inertialMat); |
323 |
|
|
324 |
|
//read dipole moment |
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 { |
416 |
|
|
417 |
|
break; |
418 |
|
|
419 |
< |
case "CubicBondType" : |
419 |
> |
case "Cubic" : |
420 |
|
if (nTokens < 4) { |
421 |
|
|
422 |
|
} else { |
430 |
|
} |
431 |
|
break; |
432 |
|
|
433 |
< |
case "QuadraticBondType" : |
433 |
> |
case "Quartic" : |
434 |
|
if (nTokens < 5) { |
435 |
|
|
436 |
|
} else { |
446 |
|
} |
447 |
|
break; |
448 |
|
|
449 |
< |
case "PolynomialBondType " : |
449 |
> |
case "Polynomial" : |
450 |
|
if (nTokens < 2 || nTokens % 2 != 0) { |
451 |
|
|
452 |
|
} else { |
499 |
|
//switch is a maintain nightmare |
500 |
|
switch(bt) { |
501 |
|
|
502 |
< |
case "HarmonicBendType" : |
502 |
> |
case "Harmonic" : |
503 |
|
|
504 |
|
if (nTokens < 1) { |
505 |
|
|
509 |
|
bendType = new HarmonicBendType(theta0, ktheta); |
510 |
|
} |
511 |
|
break; |
512 |
< |
case "GhostBendType" : |
512 |
> |
case "GhostBend" : |
513 |
|
if (nTokens < 1) { |
514 |
|
|
515 |
|
} else { |
518 |
|
} |
519 |
|
break; |
520 |
|
|
521 |
< |
case "UreyBradleyBendType" : |
521 |
> |
case "UreyBradley" : |
522 |
|
if (nTokens < 3) { |
523 |
|
|
524 |
|
} else { |
529 |
|
} |
530 |
|
break; |
531 |
|
|
532 |
< |
case "CubicBendType" : |
532 |
> |
case "Cubic" : |
533 |
|
if (nTokens < 4) { |
534 |
|
|
535 |
|
} else { |
543 |
|
} |
544 |
|
break; |
545 |
|
|
546 |
< |
case "QuadraticBendType" : |
546 |
> |
case "Quartic" : |
547 |
|
if (nTokens < 5) { |
548 |
|
|
549 |
|
} else { |
559 |
|
} |
560 |
|
break; |
561 |
|
|
562 |
< |
case "PolynomialBendType " : |
562 |
> |
case "Polynomial" : |
563 |
|
if (nTokens < 2 || nTokens % 2 != 0) { |
564 |
|
|
565 |
|
} else { |
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 |
|
|
613 |
|
|
614 |
|
switch(tt) { |
615 |
|
|
616 |
< |
case "CubicTorsionType" : |
616 |
> |
case "Cubic" : |
617 |
|
if (nTokens < 4) { |
618 |
|
|
619 |
|
} else { |
627 |
|
} |
628 |
|
break; |
629 |
|
|
630 |
< |
case "QuadraticTorsionType" : |
630 |
> |
case "Quartic" : |
631 |
|
if (nTokens < 5) { |
632 |
|
|
633 |
|
} else { |
643 |
|
} |
644 |
|
break; |
645 |
|
|
646 |
< |
case "PolynomialTorsionType " : |
646 |
> |
case "Polynomial" : |
647 |
|
if (nTokens < 2 || nTokens % 2 != 0) { |
648 |
|
|
649 |
|
} else { |
660 |
|
} |
661 |
|
|
662 |
|
break; |
663 |
< |
case "CharmmTorsionType" : |
663 |
> |
case "Charmm" : |
664 |
|
|
665 |
|
if (nTokens < 3 || nTokens % 3 != 0) { |
666 |
|
|
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 |