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 |
|
|
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 |
|
|
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; |
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 |
|
|
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); |
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 { |
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 |
|
|
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 { |
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 |
|
|
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 { |
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 |
– |
|