6 |
|
* redistribute this software in source and binary code form, provided |
7 |
|
* that the following conditions are met: |
8 |
|
* |
9 |
< |
* 1. Acknowledgement of the program authors must be made in any |
10 |
< |
* publication of scientific results based in part on use of the |
11 |
< |
* program. An acceptable form of acknowledgement is citation of |
12 |
< |
* the article in which the program was described (Matthew |
13 |
< |
* A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher |
14 |
< |
* J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented |
15 |
< |
* Parallel Simulation Engine for Molecular Dynamics," |
16 |
< |
* J. Comput. Chem. 26, pp. 252-271 (2005)) |
17 |
< |
* |
18 |
< |
* 2. Redistributions of source code must retain the above copyright |
9 |
> |
* 1. Redistributions of source code must retain the above copyright |
10 |
|
* notice, this list of conditions and the following disclaimer. |
11 |
|
* |
12 |
< |
* 3. Redistributions in binary form must reproduce the above copyright |
12 |
> |
* 2. Redistributions in binary form must reproduce the above copyright |
13 |
|
* notice, this list of conditions and the following disclaimer in the |
14 |
|
* documentation and/or other materials provided with the |
15 |
|
* distribution. |
28 |
|
* arising out of the use of or inability to use software, even if the |
29 |
|
* University of Notre Dame has been advised of the possibility of |
30 |
|
* such damages. |
31 |
+ |
* |
32 |
+ |
* SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your |
33 |
+ |
* research, please cite the appropriate papers when you publish your |
34 |
+ |
* work. Good starting points are: |
35 |
+ |
* |
36 |
+ |
* [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). |
37 |
+ |
* [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). |
38 |
+ |
* [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). |
39 |
+ |
* [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
+ |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
|
*/ |
42 |
|
|
43 |
|
/** |
48 |
|
* @version 1.0 |
49 |
|
*/ |
50 |
|
|
51 |
+ |
#include <algorithm> |
52 |
|
#include "UseTheForce/ForceField.hpp" |
53 |
|
#include "utils/simError.h" |
54 |
< |
#include "UseTheForce/DarkSide/atype_interface.h" |
55 |
< |
#include "UseTheForce/DarkSide/fForceOptions_interface.h" |
54 |
< |
#include "UseTheForce/DarkSide/switcheroo_interface.h" |
55 |
< |
namespace oopse { |
54 |
> |
#include "utils/Tuple.hpp" |
55 |
> |
#include "types/LennardJonesAdapter.hpp" |
56 |
|
|
57 |
+ |
namespace OpenMD { |
58 |
+ |
|
59 |
|
ForceField::ForceField() { |
60 |
+ |
|
61 |
|
char* tempPath; |
62 |
|
tempPath = getenv("FORCE_PARAM_PATH"); |
63 |
< |
|
63 |
> |
|
64 |
|
if (tempPath == NULL) { |
65 |
|
//convert a macro from compiler to a string in c++ |
66 |
|
STR_DEFINE(ffPath_, FRC_PATH ); |
69 |
|
} |
70 |
|
} |
71 |
|
|
72 |
< |
|
73 |
< |
ForceField::~ForceField() { |
74 |
< |
deleteAtypes(); |
75 |
< |
deleteSwitch(); |
76 |
< |
} |
77 |
< |
|
72 |
> |
/** |
73 |
> |
* getAtomType by string |
74 |
> |
* |
75 |
> |
* finds the requested atom type in this force field using the string |
76 |
> |
* name of the atom type. |
77 |
> |
*/ |
78 |
|
AtomType* ForceField::getAtomType(const std::string &at) { |
79 |
|
std::vector<std::string> keys; |
80 |
|
keys.push_back(at); |
81 |
|
return atomTypeCont_.find(keys); |
82 |
|
} |
83 |
|
|
84 |
< |
BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) { |
84 |
> |
/** |
85 |
> |
* getAtomType by ident |
86 |
> |
* |
87 |
> |
* finds the requested atom type in this force field using the |
88 |
> |
* integer ident instead of the string name of the atom type. |
89 |
> |
*/ |
90 |
> |
AtomType* ForceField::getAtomType(int ident) { |
91 |
> |
std::string at = atypeIdentToName.find(ident)->second; |
92 |
> |
return getAtomType(at); |
93 |
> |
} |
94 |
> |
|
95 |
> |
BondType* ForceField::getBondType(const std::string &at1, |
96 |
> |
const std::string &at2) { |
97 |
|
std::vector<std::string> keys; |
98 |
|
keys.push_back(at1); |
99 |
|
keys.push_back(at2); |
103 |
|
if (bondType) { |
104 |
|
return bondType; |
105 |
|
} else { |
106 |
< |
//if no exact match found, try wild card match |
107 |
< |
return bondTypeCont_.find(keys, wildCardAtomTypeName_); |
108 |
< |
} |
106 |
> |
AtomType* atype1; |
107 |
> |
AtomType* atype2; |
108 |
> |
std::vector<std::string> at1key; |
109 |
> |
at1key.push_back(at1); |
110 |
> |
atype1 = atomTypeCont_.find(at1key); |
111 |
> |
|
112 |
> |
std::vector<std::string> at2key; |
113 |
> |
at2key.push_back(at2); |
114 |
> |
atype2 = atomTypeCont_.find(at2key); |
115 |
|
|
116 |
< |
} |
116 |
> |
// query atom types for their chains of responsibility |
117 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
118 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
119 |
|
|
120 |
< |
BendType* ForceField::getBendType(const std::string &at1, const std::string &at2, |
120 |
> |
std::vector<AtomType*>::iterator i; |
121 |
> |
std::vector<AtomType*>::iterator j; |
122 |
> |
|
123 |
> |
int ii = 0; |
124 |
> |
int jj = 0; |
125 |
> |
int bondTypeScore; |
126 |
> |
|
127 |
> |
std::vector<std::pair<int, std::vector<std::string> > > foundBonds; |
128 |
> |
|
129 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
130 |
> |
jj = 0; |
131 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
132 |
> |
|
133 |
> |
bondTypeScore = ii + jj; |
134 |
> |
|
135 |
> |
std::vector<std::string> myKeys; |
136 |
> |
myKeys.push_back((*i)->getName()); |
137 |
> |
myKeys.push_back((*j)->getName()); |
138 |
> |
|
139 |
> |
BondType* bondType = bondTypeCont_.find(myKeys); |
140 |
> |
if (bondType) { |
141 |
> |
foundBonds.push_back(std::make_pair(bondTypeScore, myKeys)); |
142 |
> |
} |
143 |
> |
jj++; |
144 |
> |
} |
145 |
> |
ii++; |
146 |
> |
} |
147 |
> |
|
148 |
> |
|
149 |
> |
if (foundBonds.size() > 0) { |
150 |
> |
// sort the foundBonds by the score: |
151 |
> |
std::sort(foundBonds.begin(), foundBonds.end()); |
152 |
> |
|
153 |
> |
int bestScore = foundBonds[0].first; |
154 |
> |
std::vector<std::string> theKeys = foundBonds[0].second; |
155 |
> |
|
156 |
> |
BondType* bestType = bondTypeCont_.find(theKeys); |
157 |
> |
|
158 |
> |
return bestType; |
159 |
> |
} else { |
160 |
> |
//if no exact match found, try wild card match |
161 |
> |
return bondTypeCont_.find(keys, wildCardAtomTypeName_); |
162 |
> |
} |
163 |
> |
} |
164 |
> |
} |
165 |
> |
|
166 |
> |
BendType* ForceField::getBendType(const std::string &at1, |
167 |
> |
const std::string &at2, |
168 |
|
const std::string &at3) { |
169 |
|
std::vector<std::string> keys; |
170 |
|
keys.push_back(at1); |
176 |
|
if (bendType) { |
177 |
|
return bendType; |
178 |
|
} else { |
179 |
< |
//if no exact match found, try wild card match |
180 |
< |
return bendTypeCont_.find(keys, wildCardAtomTypeName_); |
179 |
> |
|
180 |
> |
AtomType* atype1; |
181 |
> |
AtomType* atype2; |
182 |
> |
AtomType* atype3; |
183 |
> |
std::vector<std::string> at1key; |
184 |
> |
at1key.push_back(at1); |
185 |
> |
atype1 = atomTypeCont_.find(at1key); |
186 |
> |
|
187 |
> |
std::vector<std::string> at2key; |
188 |
> |
at2key.push_back(at2); |
189 |
> |
atype2 = atomTypeCont_.find(at2key); |
190 |
> |
|
191 |
> |
std::vector<std::string> at3key; |
192 |
> |
at3key.push_back(at3); |
193 |
> |
atype3 = atomTypeCont_.find(at3key); |
194 |
> |
|
195 |
> |
// query atom types for their chains of responsibility |
196 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
197 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
198 |
> |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
199 |
> |
|
200 |
> |
std::vector<AtomType*>::iterator i; |
201 |
> |
std::vector<AtomType*>::iterator j; |
202 |
> |
std::vector<AtomType*>::iterator k; |
203 |
> |
|
204 |
> |
int ii = 0; |
205 |
> |
int jj = 0; |
206 |
> |
int kk = 0; |
207 |
> |
int IKscore; |
208 |
> |
|
209 |
> |
std::vector<tuple3<int, int, std::vector<std::string> > > foundBends; |
210 |
> |
|
211 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
212 |
> |
ii = 0; |
213 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
214 |
> |
kk = 0; |
215 |
> |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
216 |
> |
|
217 |
> |
IKscore = ii + kk; |
218 |
> |
|
219 |
> |
std::vector<std::string> myKeys; |
220 |
> |
myKeys.push_back((*i)->getName()); |
221 |
> |
myKeys.push_back((*j)->getName()); |
222 |
> |
myKeys.push_back((*k)->getName()); |
223 |
> |
|
224 |
> |
BendType* bendType = bendTypeCont_.find(myKeys); |
225 |
> |
if (bendType) { |
226 |
> |
foundBends.push_back( make_tuple3(jj, IKscore, myKeys) ); |
227 |
> |
} |
228 |
> |
kk++; |
229 |
> |
} |
230 |
> |
ii++; |
231 |
> |
} |
232 |
> |
jj++; |
233 |
> |
} |
234 |
> |
|
235 |
> |
if (foundBends.size() > 0) { |
236 |
> |
std::sort(foundBends.begin(), foundBends.end()); |
237 |
> |
int jscore = foundBends[0].first; |
238 |
> |
int ikscore = foundBends[0].second; |
239 |
> |
std::vector<std::string> theKeys = foundBends[0].third; |
240 |
> |
|
241 |
> |
BendType* bestType = bendTypeCont_.find(theKeys); |
242 |
> |
return bestType; |
243 |
> |
} else { |
244 |
> |
//if no exact match found, try wild card match |
245 |
> |
return bendTypeCont_.find(keys, wildCardAtomTypeName_); |
246 |
> |
} |
247 |
|
} |
248 |
|
} |
249 |
|
|
250 |
< |
TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2, |
251 |
< |
const std::string &at3, const std::string &at4) { |
250 |
> |
TorsionType* ForceField::getTorsionType(const std::string &at1, |
251 |
> |
const std::string &at2, |
252 |
> |
const std::string &at3, |
253 |
> |
const std::string &at4) { |
254 |
|
std::vector<std::string> keys; |
255 |
|
keys.push_back(at1); |
256 |
|
keys.push_back(at2); |
257 |
|
keys.push_back(at3); |
258 |
|
keys.push_back(at4); |
259 |
|
|
260 |
+ |
|
261 |
+ |
//try exact match first |
262 |
|
TorsionType* torsionType = torsionTypeCont_.find(keys); |
263 |
|
if (torsionType) { |
264 |
|
return torsionType; |
265 |
|
} else { |
126 |
– |
//if no exact match found, try wild card match |
127 |
– |
return torsionTypeCont_.find(keys, wildCardAtomTypeName_); |
128 |
– |
} |
129 |
– |
|
130 |
– |
return torsionTypeCont_.find(keys, wildCardAtomTypeName_); |
266 |
|
|
267 |
+ |
AtomType* atype1; |
268 |
+ |
AtomType* atype2; |
269 |
+ |
AtomType* atype3; |
270 |
+ |
AtomType* atype4; |
271 |
+ |
std::vector<std::string> at1key; |
272 |
+ |
at1key.push_back(at1); |
273 |
+ |
atype1 = atomTypeCont_.find(at1key); |
274 |
+ |
|
275 |
+ |
std::vector<std::string> at2key; |
276 |
+ |
at2key.push_back(at2); |
277 |
+ |
atype2 = atomTypeCont_.find(at2key); |
278 |
+ |
|
279 |
+ |
std::vector<std::string> at3key; |
280 |
+ |
at3key.push_back(at3); |
281 |
+ |
atype3 = atomTypeCont_.find(at3key); |
282 |
+ |
|
283 |
+ |
std::vector<std::string> at4key; |
284 |
+ |
at4key.push_back(at4); |
285 |
+ |
atype4 = atomTypeCont_.find(at4key); |
286 |
+ |
|
287 |
+ |
// query atom types for their chains of responsibility |
288 |
+ |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
289 |
+ |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
290 |
+ |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
291 |
+ |
std::vector<AtomType*> at4Chain = atype4->allYourBase(); |
292 |
+ |
|
293 |
+ |
std::vector<AtomType*>::iterator i; |
294 |
+ |
std::vector<AtomType*>::iterator j; |
295 |
+ |
std::vector<AtomType*>::iterator k; |
296 |
+ |
std::vector<AtomType*>::iterator l; |
297 |
+ |
|
298 |
+ |
int ii = 0; |
299 |
+ |
int jj = 0; |
300 |
+ |
int kk = 0; |
301 |
+ |
int ll = 0; |
302 |
+ |
int ILscore; |
303 |
+ |
int JKscore; |
304 |
+ |
|
305 |
+ |
std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions; |
306 |
+ |
|
307 |
+ |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
308 |
+ |
kk = 0; |
309 |
+ |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
310 |
+ |
ii = 0; |
311 |
+ |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
312 |
+ |
ll = 0; |
313 |
+ |
for (l = at4Chain.begin(); l != at4Chain.end(); l++) { |
314 |
+ |
|
315 |
+ |
ILscore = ii + ll; |
316 |
+ |
JKscore = jj + kk; |
317 |
+ |
|
318 |
+ |
std::vector<std::string> myKeys; |
319 |
+ |
myKeys.push_back((*i)->getName()); |
320 |
+ |
myKeys.push_back((*j)->getName()); |
321 |
+ |
myKeys.push_back((*k)->getName()); |
322 |
+ |
myKeys.push_back((*l)->getName()); |
323 |
+ |
|
324 |
+ |
TorsionType* torsionType = torsionTypeCont_.find(myKeys); |
325 |
+ |
if (torsionType) { |
326 |
+ |
foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) ); |
327 |
+ |
} |
328 |
+ |
ll++; |
329 |
+ |
} |
330 |
+ |
ii++; |
331 |
+ |
} |
332 |
+ |
kk++; |
333 |
+ |
} |
334 |
+ |
jj++; |
335 |
+ |
} |
336 |
+ |
|
337 |
+ |
if (foundTorsions.size() > 0) { |
338 |
+ |
std::sort(foundTorsions.begin(), foundTorsions.end()); |
339 |
+ |
int jkscore = foundTorsions[0].first; |
340 |
+ |
int ilscore = foundTorsions[0].second; |
341 |
+ |
std::vector<std::string> theKeys = foundTorsions[0].third; |
342 |
+ |
|
343 |
+ |
TorsionType* bestType = torsionTypeCont_.find(theKeys); |
344 |
+ |
return bestType; |
345 |
+ |
} else { |
346 |
+ |
//if no exact match found, try wild card match |
347 |
+ |
return torsionTypeCont_.find(keys, wildCardAtomTypeName_); |
348 |
+ |
} |
349 |
+ |
} |
350 |
|
} |
351 |
|
|
352 |
< |
BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){ |
352 |
> |
InversionType* ForceField::getInversionType(const std::string &at1, |
353 |
> |
const std::string &at2, |
354 |
> |
const std::string &at3, |
355 |
> |
const std::string &at4) { |
356 |
> |
std::vector<std::string> keys; |
357 |
> |
keys.push_back(at1); |
358 |
> |
keys.push_back(at2); |
359 |
> |
keys.push_back(at3); |
360 |
> |
keys.push_back(at4); |
361 |
> |
|
362 |
> |
//try exact match first |
363 |
> |
InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); |
364 |
> |
if (inversionType) { |
365 |
> |
return inversionType; |
366 |
> |
} else { |
367 |
> |
|
368 |
> |
AtomType* atype1; |
369 |
> |
AtomType* atype2; |
370 |
> |
AtomType* atype3; |
371 |
> |
AtomType* atype4; |
372 |
> |
std::vector<std::string> at1key; |
373 |
> |
at1key.push_back(at1); |
374 |
> |
atype1 = atomTypeCont_.find(at1key); |
375 |
> |
|
376 |
> |
std::vector<std::string> at2key; |
377 |
> |
at2key.push_back(at2); |
378 |
> |
atype2 = atomTypeCont_.find(at2key); |
379 |
> |
|
380 |
> |
std::vector<std::string> at3key; |
381 |
> |
at3key.push_back(at3); |
382 |
> |
atype3 = atomTypeCont_.find(at3key); |
383 |
> |
|
384 |
> |
std::vector<std::string> at4key; |
385 |
> |
at4key.push_back(at4); |
386 |
> |
atype4 = atomTypeCont_.find(at4key); |
387 |
> |
|
388 |
> |
// query atom types for their chains of responsibility |
389 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
390 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
391 |
> |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
392 |
> |
std::vector<AtomType*> at4Chain = atype4->allYourBase(); |
393 |
> |
|
394 |
> |
std::vector<AtomType*>::iterator i; |
395 |
> |
std::vector<AtomType*>::iterator j; |
396 |
> |
std::vector<AtomType*>::iterator k; |
397 |
> |
std::vector<AtomType*>::iterator l; |
398 |
> |
|
399 |
> |
int ii = 0; |
400 |
> |
int jj = 0; |
401 |
> |
int kk = 0; |
402 |
> |
int ll = 0; |
403 |
> |
int Iscore; |
404 |
> |
int JKLscore; |
405 |
> |
|
406 |
> |
std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions; |
407 |
> |
|
408 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
409 |
> |
kk = 0; |
410 |
> |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
411 |
> |
ii = 0; |
412 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
413 |
> |
ll = 0; |
414 |
> |
for (l = at4Chain.begin(); l != at4Chain.end(); l++) { |
415 |
> |
|
416 |
> |
Iscore = ii; |
417 |
> |
JKLscore = jj + kk + ll; |
418 |
> |
|
419 |
> |
std::vector<std::string> myKeys; |
420 |
> |
myKeys.push_back((*i)->getName()); |
421 |
> |
myKeys.push_back((*j)->getName()); |
422 |
> |
myKeys.push_back((*k)->getName()); |
423 |
> |
myKeys.push_back((*l)->getName()); |
424 |
> |
|
425 |
> |
InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); |
426 |
> |
if (inversionType) { |
427 |
> |
foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); |
428 |
> |
} |
429 |
> |
ll++; |
430 |
> |
} |
431 |
> |
ii++; |
432 |
> |
} |
433 |
> |
kk++; |
434 |
> |
} |
435 |
> |
jj++; |
436 |
> |
} |
437 |
> |
|
438 |
> |
if (foundInversions.size() > 0) { |
439 |
> |
std::sort(foundInversions.begin(), foundInversions.end()); |
440 |
> |
int iscore = foundInversions[0].first; |
441 |
> |
int jklscore = foundInversions[0].second; |
442 |
> |
std::vector<std::string> theKeys = foundInversions[0].third; |
443 |
> |
|
444 |
> |
InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); |
445 |
> |
return bestType; |
446 |
> |
} else { |
447 |
> |
//if no exact match found, try wild card match |
448 |
> |
return inversionTypeCont_.find(keys, wildCardAtomTypeName_); |
449 |
> |
} |
450 |
> |
} |
451 |
> |
} |
452 |
> |
|
453 |
> |
NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { |
454 |
> |
|
455 |
|
std::vector<std::string> keys; |
456 |
|
keys.push_back(at1); |
457 |
|
keys.push_back(at2); |
458 |
+ |
|
459 |
+ |
//try exact match first |
460 |
+ |
NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); |
461 |
+ |
if (nbiType) { |
462 |
+ |
return nbiType; |
463 |
+ |
} else { |
464 |
+ |
AtomType* atype1; |
465 |
+ |
AtomType* atype2; |
466 |
+ |
std::vector<std::string> at1key; |
467 |
+ |
at1key.push_back(at1); |
468 |
+ |
atype1 = atomTypeCont_.find(at1key); |
469 |
+ |
|
470 |
+ |
std::vector<std::string> at2key; |
471 |
+ |
at2key.push_back(at2); |
472 |
+ |
atype2 = atomTypeCont_.find(at2key); |
473 |
+ |
|
474 |
+ |
// query atom types for their chains of responsibility |
475 |
+ |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
476 |
+ |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
477 |
+ |
|
478 |
+ |
std::vector<AtomType*>::iterator i; |
479 |
+ |
std::vector<AtomType*>::iterator j; |
480 |
+ |
|
481 |
+ |
int ii = 0; |
482 |
+ |
int jj = 0; |
483 |
+ |
int nbiTypeScore; |
484 |
+ |
|
485 |
+ |
std::vector<std::pair<int, std::vector<std::string> > > foundNBI; |
486 |
+ |
|
487 |
+ |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
488 |
+ |
jj = 0; |
489 |
+ |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
490 |
+ |
|
491 |
+ |
nbiTypeScore = ii + jj; |
492 |
+ |
|
493 |
+ |
std::vector<std::string> myKeys; |
494 |
+ |
myKeys.push_back((*i)->getName()); |
495 |
+ |
myKeys.push_back((*j)->getName()); |
496 |
+ |
|
497 |
+ |
NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys); |
498 |
+ |
if (nbiType) { |
499 |
+ |
foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys)); |
500 |
+ |
} |
501 |
+ |
jj++; |
502 |
+ |
} |
503 |
+ |
ii++; |
504 |
+ |
} |
505 |
+ |
|
506 |
+ |
|
507 |
+ |
if (foundNBI.size() > 0) { |
508 |
+ |
// sort the foundNBI by the score: |
509 |
+ |
std::sort(foundNBI.begin(), foundNBI.end()); |
510 |
+ |
|
511 |
+ |
int bestScore = foundNBI[0].first; |
512 |
+ |
std::vector<std::string> theKeys = foundNBI[0].second; |
513 |
+ |
|
514 |
+ |
NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys); |
515 |
+ |
return bestType; |
516 |
+ |
} else { |
517 |
+ |
//if no exact match found, try wild card match |
518 |
+ |
return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); |
519 |
+ |
} |
520 |
+ |
} |
521 |
+ |
} |
522 |
+ |
|
523 |
+ |
BondType* ForceField::getExactBondType(const std::string &at1, |
524 |
+ |
const std::string &at2){ |
525 |
+ |
std::vector<std::string> keys; |
526 |
+ |
keys.push_back(at1); |
527 |
+ |
keys.push_back(at2); |
528 |
|
return bondTypeCont_.find(keys); |
529 |
|
} |
530 |
< |
|
531 |
< |
BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2, |
530 |
> |
|
531 |
> |
BendType* ForceField::getExactBendType(const std::string &at1, |
532 |
> |
const std::string &at2, |
533 |
|
const std::string &at3){ |
534 |
|
std::vector<std::string> keys; |
535 |
|
keys.push_back(at1); |
537 |
|
keys.push_back(at3); |
538 |
|
return bendTypeCont_.find(keys); |
539 |
|
} |
540 |
< |
|
541 |
< |
TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2, |
542 |
< |
const std::string &at3, const std::string &at4){ |
540 |
> |
|
541 |
> |
TorsionType* ForceField::getExactTorsionType(const std::string &at1, |
542 |
> |
const std::string &at2, |
543 |
> |
const std::string &at3, |
544 |
> |
const std::string &at4){ |
545 |
|
std::vector<std::string> keys; |
546 |
|
keys.push_back(at1); |
547 |
|
keys.push_back(at2); |
549 |
|
keys.push_back(at4); |
550 |
|
return torsionTypeCont_.find(keys); |
551 |
|
} |
552 |
+ |
|
553 |
+ |
InversionType* ForceField::getExactInversionType(const std::string &at1, |
554 |
+ |
const std::string &at2, |
555 |
+ |
const std::string &at3, |
556 |
+ |
const std::string &at4){ |
557 |
+ |
std::vector<std::string> keys; |
558 |
+ |
keys.push_back(at1); |
559 |
+ |
keys.push_back(at2); |
560 |
+ |
keys.push_back(at3); |
561 |
+ |
keys.push_back(at4); |
562 |
+ |
return inversionTypeCont_.find(keys); |
563 |
+ |
} |
564 |
+ |
|
565 |
+ |
NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ |
566 |
+ |
std::vector<std::string> keys; |
567 |
+ |
keys.push_back(at1); |
568 |
+ |
keys.push_back(at2); |
569 |
+ |
return nonBondedInteractionTypeCont_.find(keys); |
570 |
+ |
} |
571 |
+ |
|
572 |
+ |
|
573 |
|
bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { |
574 |
|
std::vector<std::string> keys; |
575 |
|
keys.push_back(at); |
576 |
+ |
atypeIdentToName[atomType->getIdent()] = at; |
577 |
|
return atomTypeCont_.add(keys, atomType); |
578 |
|
} |
579 |
|
|
580 |
< |
bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) { |
580 |
> |
bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { |
581 |
|
std::vector<std::string> keys; |
582 |
+ |
keys.push_back(at); |
583 |
+ |
atypeIdentToName[atomType->getIdent()] = at; |
584 |
+ |
return atomTypeCont_.replace(keys, atomType); |
585 |
+ |
} |
586 |
+ |
|
587 |
+ |
bool ForceField::addBondType(const std::string &at1, const std::string &at2, |
588 |
+ |
BondType* bondType) { |
589 |
+ |
std::vector<std::string> keys; |
590 |
|
keys.push_back(at1); |
591 |
|
keys.push_back(at2); |
592 |
< |
return bondTypeCont_.add(keys, bondType); |
170 |
< |
|
592 |
> |
return bondTypeCont_.add(keys, bondType); |
593 |
|
} |
594 |
< |
|
594 |
> |
|
595 |
|
bool ForceField::addBendType(const std::string &at1, const std::string &at2, |
596 |
|
const std::string &at3, BendType* bendType) { |
597 |
|
std::vector<std::string> keys; |
600 |
|
keys.push_back(at3); |
601 |
|
return bendTypeCont_.add(keys, bendType); |
602 |
|
} |
603 |
< |
|
604 |
< |
bool ForceField::addTorsionType(const std::string &at1, const std::string &at2, |
605 |
< |
const std::string &at3, const std::string &at4, TorsionType* torsionType) { |
603 |
> |
|
604 |
> |
bool ForceField::addTorsionType(const std::string &at1, |
605 |
> |
const std::string &at2, |
606 |
> |
const std::string &at3, |
607 |
> |
const std::string &at4, |
608 |
> |
TorsionType* torsionType) { |
609 |
|
std::vector<std::string> keys; |
610 |
|
keys.push_back(at1); |
611 |
|
keys.push_back(at2); |
614 |
|
return torsionTypeCont_.add(keys, torsionType); |
615 |
|
} |
616 |
|
|
617 |
+ |
bool ForceField::addInversionType(const std::string &at1, |
618 |
+ |
const std::string &at2, |
619 |
+ |
const std::string &at3, |
620 |
+ |
const std::string &at4, |
621 |
+ |
InversionType* inversionType) { |
622 |
+ |
std::vector<std::string> keys; |
623 |
+ |
keys.push_back(at1); |
624 |
+ |
keys.push_back(at2); |
625 |
+ |
keys.push_back(at3); |
626 |
+ |
keys.push_back(at4); |
627 |
+ |
return inversionTypeCont_.add(keys, inversionType); |
628 |
+ |
} |
629 |
+ |
|
630 |
+ |
bool ForceField::addNonBondedInteractionType(const std::string &at1, |
631 |
+ |
const std::string &at2, |
632 |
+ |
NonBondedInteractionType* nbiType) { |
633 |
+ |
std::vector<std::string> keys; |
634 |
+ |
keys.push_back(at1); |
635 |
+ |
keys.push_back(at2); |
636 |
+ |
return nonBondedInteractionTypeCont_.add(keys, nbiType); |
637 |
+ |
} |
638 |
+ |
|
639 |
|
RealType ForceField::getRcutFromAtomType(AtomType* at) { |
193 |
– |
/**@todo */ |
194 |
– |
GenericData* data; |
640 |
|
RealType rcut = 0.0; |
641 |
< |
|
642 |
< |
if (at->isLennardJones()) { |
643 |
< |
data = at->getPropertyByName("LennardJones"); |
644 |
< |
if (data != NULL) { |
200 |
< |
LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); |
201 |
< |
|
202 |
< |
if (ljData != NULL) { |
203 |
< |
LJParam ljParam = ljData->getData(); |
204 |
< |
|
205 |
< |
//by default use 2.5*sigma as cutoff radius |
206 |
< |
rcut = 2.5 * ljParam.sigma; |
207 |
< |
|
208 |
< |
} else { |
209 |
< |
sprintf( painCave.errMsg, |
210 |
< |
"Can not cast GenericData to LJParam\n"); |
211 |
< |
painCave.severity = OOPSE_ERROR; |
212 |
< |
painCave.isFatal = 1; |
213 |
< |
simError(); |
214 |
< |
} |
215 |
< |
} else { |
216 |
< |
sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); |
217 |
< |
painCave.severity = OOPSE_ERROR; |
218 |
< |
painCave.isFatal = 1; |
219 |
< |
simError(); |
220 |
< |
} |
641 |
> |
|
642 |
> |
LennardJonesAdapter lja = LennardJonesAdapter(at); |
643 |
> |
if (lja.isLennardJones()) { |
644 |
> |
rcut = 2.5 * lja.getSigma(); |
645 |
|
} |
222 |
– |
|
646 |
|
return rcut; |
647 |
|
} |
648 |
+ |
|
649 |
|
|
226 |
– |
|
650 |
|
ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { |
651 |
|
std::string forceFieldFilename(filename); |
652 |
|
ifstrstream* ffStream = new ifstrstream(); |
668 |
|
"\tHave you tried setting the FORCE_PARAM_PATH environment " |
669 |
|
"variable?\n", |
670 |
|
forceFieldFilename.c_str() ); |
671 |
< |
painCave.severity = OOPSE_ERROR; |
671 |
> |
painCave.severity = OPENMD_ERROR; |
672 |
|
painCave.isFatal = 1; |
673 |
|
simError(); |
674 |
|
} |
675 |
|
} |
253 |
– |
|
676 |
|
return ffStream; |
255 |
– |
|
677 |
|
} |
678 |
|
|
679 |
< |
void ForceField::setFortranForceOptions(){ |
259 |
< |
ForceOptions theseFortranOptions; |
260 |
< |
forceFieldOptions_.makeFortranOptions(theseFortranOptions); |
261 |
< |
setfForceOptions(&theseFortranOptions); |
262 |
< |
} |
263 |
< |
} //end namespace oopse |
679 |
> |
} //end namespace OpenMD |