1 |
+ |
/* |
2 |
+ |
* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved. |
3 |
+ |
* |
4 |
+ |
* The University of Notre Dame grants you ("Licensee") a |
5 |
+ |
* non-exclusive, royalty free, license to use, modify and |
6 |
+ |
* redistribute this software in source and binary code form, provided |
7 |
+ |
* that the following conditions are met: |
8 |
+ |
* |
9 |
+ |
* 1. Redistributions of source code must retain the above copyright |
10 |
+ |
* notice, this list of conditions and the following disclaimer. |
11 |
+ |
* |
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. |
16 |
+ |
* |
17 |
+ |
* This software is provided "AS IS," without a warranty of any |
18 |
+ |
* kind. All express or implied conditions, representations and |
19 |
+ |
* warranties, including any implied warranty of merchantability, |
20 |
+ |
* fitness for a particular purpose or non-infringement, are hereby |
21 |
+ |
* excluded. The University of Notre Dame and its licensors shall not |
22 |
+ |
* be liable for any damages suffered by licensee as a result of |
23 |
+ |
* using, modifying or distributing the software or its |
24 |
+ |
* derivatives. In no event will the University of Notre Dame or its |
25 |
+ |
* licensors be liable for any lost revenue, profit or data, or for |
26 |
+ |
* direct, indirect, special, consequential, incidental or punitive |
27 |
+ |
* damages, however caused and regardless of the theory of liability, |
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] Vardeman & Gezelter, in progress (2009). |
40 |
+ |
*/ |
41 |
+ |
|
42 |
+ |
/** |
43 |
+ |
* @file ForceField.cpp |
44 |
+ |
* @author tlin |
45 |
+ |
* @date 11/04/2004 |
46 |
+ |
* @time 22:51am |
47 |
+ |
* @version 1.0 |
48 |
+ |
*/ |
49 |
+ |
|
50 |
+ |
#include <algorithm> |
51 |
|
#include "UseTheForce/ForceField.hpp" |
52 |
+ |
#include "utils/simError.h" |
53 |
+ |
#include "utils/Tuple.hpp" |
54 |
+ |
#include "UseTheForce/DarkSide/atype_interface.h" |
55 |
+ |
namespace OpenMD { |
56 |
|
|
57 |
< |
AtomType* ForceField::getMatchingAtomType(const string &at) { |
57 |
> |
ForceField::ForceField() { |
58 |
|
|
59 |
< |
map<string, AtomType*>::iterator iter; |
60 |
< |
|
61 |
< |
iter = atomTypeMap.find(at); |
62 |
< |
if (iter != atomTypeMap.end()) { |
63 |
< |
return iter->second; |
64 |
< |
} else { |
65 |
< |
return NULL; |
59 |
> |
char* tempPath; |
60 |
> |
tempPath = getenv("FORCE_PARAM_PATH"); |
61 |
> |
|
62 |
> |
if (tempPath == NULL) { |
63 |
> |
//convert a macro from compiler to a string in c++ |
64 |
> |
STR_DEFINE(ffPath_, FRC_PATH ); |
65 |
> |
} else { |
66 |
> |
ffPath_ = tempPath; |
67 |
> |
} |
68 |
|
} |
13 |
– |
} |
69 |
|
|
15 |
– |
BondType* ForceField::getMatchingBondType(const string &at1, |
16 |
– |
const string &at2) { |
70 |
|
|
71 |
< |
map<pair<string,string>, BondType*>::iterator iter; |
72 |
< |
vector<BondType*> foundTypes; |
71 |
> |
ForceField::~ForceField() { |
72 |
> |
deleteAtypes(); |
73 |
> |
} |
74 |
|
|
75 |
< |
iter = bondTypeMap.find(pair<at1, at2>); |
76 |
< |
if (iter != bondTypeMap.end()) { |
77 |
< |
// exact match, so just return it |
78 |
< |
return iter->second; |
79 |
< |
} |
75 |
> |
AtomType* ForceField::getAtomType(const std::string &at) { |
76 |
> |
std::vector<std::string> keys; |
77 |
> |
keys.push_back(at); |
78 |
> |
return atomTypeCont_.find(keys); |
79 |
> |
} |
80 |
|
|
81 |
< |
iter = bondTypeMap.find(pair<at2, at1>); |
82 |
< |
if (iter != bondTypeMap.end()) { |
83 |
< |
// exact match in reverse order, so just return it |
84 |
< |
return iter->second; |
85 |
< |
} |
81 |
> |
BondType* ForceField::getBondType(const std::string &at1, |
82 |
> |
const std::string &at2) { |
83 |
> |
std::vector<std::string> keys; |
84 |
> |
keys.push_back(at1); |
85 |
> |
keys.push_back(at2); |
86 |
|
|
87 |
< |
iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>); |
88 |
< |
if (iter != bondTypeMap.end()) { |
89 |
< |
foundTypes.push_back(iter->second); |
90 |
< |
} |
87 |
> |
//try exact match first |
88 |
> |
BondType* bondType = bondTypeCont_.find(keys); |
89 |
> |
if (bondType) { |
90 |
> |
return bondType; |
91 |
> |
} else { |
92 |
> |
AtomType* atype1; |
93 |
> |
AtomType* atype2; |
94 |
> |
std::vector<std::string> at1key; |
95 |
> |
at1key.push_back(at1); |
96 |
> |
atype1 = atomTypeCont_.find(at1key); |
97 |
> |
|
98 |
> |
std::vector<std::string> at2key; |
99 |
> |
at2key.push_back(at2); |
100 |
> |
atype2 = atomTypeCont_.find(at2key); |
101 |
|
|
102 |
< |
iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>); |
103 |
< |
if (iter != bondTypeMap.end()) { |
104 |
< |
foundTypes.push_back(iter->second); |
41 |
< |
} |
102 |
> |
// query atom types for their chains of responsibility |
103 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
104 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
105 |
|
|
106 |
< |
iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>); |
107 |
< |
if (iter != bondTypeMap.end()) { |
45 |
< |
foundTypes.push_back(iter->second); |
46 |
< |
} |
106 |
> |
std::vector<AtomType*>::iterator i; |
107 |
> |
std::vector<AtomType*>::iterator j; |
108 |
|
|
109 |
< |
iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>); |
110 |
< |
if (iter != bondTypeMap.end()) { |
111 |
< |
foundTypes.push_back(iter->second); |
109 |
> |
int ii = 0; |
110 |
> |
int jj = 0; |
111 |
> |
int bondTypeScore; |
112 |
> |
|
113 |
> |
std::vector<std::pair<int, std::vector<std::string> > > foundBonds; |
114 |
> |
|
115 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
116 |
> |
jj = 0; |
117 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
118 |
> |
|
119 |
> |
bondTypeScore = ii + jj; |
120 |
> |
|
121 |
> |
std::vector<std::string> myKeys; |
122 |
> |
myKeys.push_back((*i)->getName()); |
123 |
> |
myKeys.push_back((*j)->getName()); |
124 |
> |
|
125 |
> |
BondType* bondType = bondTypeCont_.find(myKeys); |
126 |
> |
if (bondType) { |
127 |
> |
foundBonds.push_back(std::make_pair(bondTypeScore, myKeys)); |
128 |
> |
} |
129 |
> |
jj++; |
130 |
> |
} |
131 |
> |
ii++; |
132 |
> |
} |
133 |
> |
|
134 |
> |
|
135 |
> |
if (foundBonds.size() > 0) { |
136 |
> |
// sort the foundBonds by the score: |
137 |
> |
std::sort(foundBonds.begin(), foundBonds.end()); |
138 |
> |
|
139 |
> |
int bestScore = foundBonds[0].first; |
140 |
> |
std::vector<std::string> theKeys = foundBonds[0].second; |
141 |
> |
|
142 |
> |
BondType* bestType = bondTypeCont_.find(theKeys); |
143 |
> |
|
144 |
> |
return bestType; |
145 |
> |
} else { |
146 |
> |
//if no exact match found, try wild card match |
147 |
> |
return bondTypeCont_.find(keys, wildCardAtomTypeName_); |
148 |
> |
} |
149 |
> |
} |
150 |
|
} |
151 |
|
|
152 |
< |
if (foundTypes.empty()) { |
153 |
< |
return NULL; |
154 |
< |
} else { |
155 |
< |
|
152 |
> |
BendType* ForceField::getBendType(const std::string &at1, |
153 |
> |
const std::string &at2, |
154 |
> |
const std::string &at3) { |
155 |
> |
std::vector<std::string> keys; |
156 |
> |
keys.push_back(at1); |
157 |
> |
keys.push_back(at2); |
158 |
> |
keys.push_back(at3); |
159 |
|
|
160 |
< |
|
160 |
> |
//try exact match first |
161 |
> |
BendType* bendType = bendTypeCont_.find(keys); |
162 |
> |
if (bendType) { |
163 |
> |
return bendType; |
164 |
> |
} else { |
165 |
|
|
166 |
+ |
AtomType* atype1; |
167 |
+ |
AtomType* atype2; |
168 |
+ |
AtomType* atype3; |
169 |
+ |
std::vector<std::string> at1key; |
170 |
+ |
at1key.push_back(at1); |
171 |
+ |
atype1 = atomTypeCont_.find(at1key); |
172 |
+ |
|
173 |
+ |
std::vector<std::string> at2key; |
174 |
+ |
at2key.push_back(at2); |
175 |
+ |
atype2 = atomTypeCont_.find(at2key); |
176 |
|
|
177 |
+ |
std::vector<std::string> at3key; |
178 |
+ |
at3key.push_back(at3); |
179 |
+ |
atype3 = atomTypeCont_.find(at3key); |
180 |
|
|
181 |
< |
|
181 |
> |
// query atom types for their chains of responsibility |
182 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
183 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
184 |
> |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
185 |
|
|
186 |
+ |
std::vector<AtomType*>::iterator i; |
187 |
+ |
std::vector<AtomType*>::iterator j; |
188 |
+ |
std::vector<AtomType*>::iterator k; |
189 |
|
|
190 |
< |
BendType* ForceField::getMatchingBendType(const string &at1, const string &at2, |
191 |
< |
const string &at3); |
192 |
< |
TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2, |
193 |
< |
const string &at3, const string &at4); |
190 |
> |
int ii = 0; |
191 |
> |
int jj = 0; |
192 |
> |
int kk = 0; |
193 |
> |
int IKscore; |
194 |
|
|
195 |
< |
double ForceField::getRcutForAtomType(AtomType* at); |
195 |
> |
std::vector<tuple3<int, int, std::vector<std::string> > > foundBends; |
196 |
|
|
197 |
+ |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
198 |
+ |
ii = 0; |
199 |
+ |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
200 |
+ |
kk = 0; |
201 |
+ |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
202 |
+ |
|
203 |
+ |
IKscore = ii + kk; |
204 |
|
|
205 |
< |
vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) { |
206 |
< |
|
207 |
< |
vector<vector<string> > results; |
205 |
> |
std::vector<std::string> myKeys; |
206 |
> |
myKeys.push_back((*i)->getName()); |
207 |
> |
myKeys.push_back((*j)->getName()); |
208 |
> |
myKeys.push_back((*k)->getName()); |
209 |
|
|
210 |
< |
|
210 |
> |
BendType* bendType = bendTypeCont_.find(myKeys); |
211 |
> |
if (bendType) { |
212 |
> |
foundBends.push_back( make_tuple3(jj, IKscore, myKeys) ); |
213 |
> |
} |
214 |
> |
kk++; |
215 |
> |
} |
216 |
> |
ii++; |
217 |
> |
} |
218 |
> |
jj++; |
219 |
> |
} |
220 |
> |
|
221 |
> |
if (foundBends.size() > 0) { |
222 |
> |
std::sort(foundBends.begin(), foundBends.end()); |
223 |
> |
int jscore = foundBends[0].first; |
224 |
> |
int ikscore = foundBends[0].second; |
225 |
> |
std::vector<std::string> theKeys = foundBends[0].third; |
226 |
> |
|
227 |
> |
BendType* bestType = bendTypeCont_.find(theKeys); |
228 |
> |
return bestType; |
229 |
> |
} else { |
230 |
> |
//if no exact match found, try wild card match |
231 |
> |
return bendTypeCont_.find(keys, wildCardAtomTypeName_); |
232 |
> |
} |
233 |
> |
} |
234 |
> |
} |
235 |
|
|
236 |
+ |
TorsionType* ForceField::getTorsionType(const std::string &at1, |
237 |
+ |
const std::string &at2, |
238 |
+ |
const std::string &at3, |
239 |
+ |
const std::string &at4) { |
240 |
+ |
std::vector<std::string> keys; |
241 |
+ |
keys.push_back(at1); |
242 |
+ |
keys.push_back(at2); |
243 |
+ |
keys.push_back(at3); |
244 |
+ |
keys.push_back(at4); |
245 |
|
|
80 |
– |
vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) { |
81 |
– |
|
82 |
– |
int nStrings; |
83 |
– |
vector<string> oneResult; |
84 |
– |
vector<vector<string> > allResults; |
246 |
|
|
247 |
< |
nStrings = myAts.size(); |
247 |
> |
//try exact match first |
248 |
> |
TorsionType* torsionType = torsionTypeCont_.find(keys); |
249 |
> |
if (torsionType) { |
250 |
> |
return torsionType; |
251 |
> |
} else { |
252 |
|
|
253 |
< |
if (nStrings == 1) { |
254 |
< |
oneResult.push_back(wildcardCharacter); |
255 |
< |
allResults.push_back(oneResult); |
256 |
< |
return allResults; |
257 |
< |
} else { |
258 |
< |
|
259 |
< |
for (i=0; i < nStrings; i++) { |
260 |
< |
oneResult = myAts; |
261 |
< |
replace(oneResult.begin(), oneResult.end(), |
253 |
> |
AtomType* atype1; |
254 |
> |
AtomType* atype2; |
255 |
> |
AtomType* atype3; |
256 |
> |
AtomType* atype4; |
257 |
> |
std::vector<std::string> at1key; |
258 |
> |
at1key.push_back(at1); |
259 |
> |
atype1 = atomTypeCont_.find(at1key); |
260 |
> |
|
261 |
> |
std::vector<std::string> at2key; |
262 |
> |
at2key.push_back(at2); |
263 |
> |
atype2 = atomTypeCont_.find(at2key); |
264 |
> |
|
265 |
> |
std::vector<std::string> at3key; |
266 |
> |
at3key.push_back(at3); |
267 |
> |
atype3 = atomTypeCont_.find(at3key); |
268 |
> |
|
269 |
> |
std::vector<std::string> at4key; |
270 |
> |
at4key.push_back(at4); |
271 |
> |
atype4 = atomTypeCont_.find(at4key); |
272 |
> |
|
273 |
> |
// query atom types for their chains of responsibility |
274 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
275 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
276 |
> |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
277 |
> |
std::vector<AtomType*> at4Chain = atype4->allYourBase(); |
278 |
> |
|
279 |
> |
std::vector<AtomType*>::iterator i; |
280 |
> |
std::vector<AtomType*>::iterator j; |
281 |
> |
std::vector<AtomType*>::iterator k; |
282 |
> |
std::vector<AtomType*>::iterator l; |
283 |
> |
|
284 |
> |
int ii = 0; |
285 |
> |
int jj = 0; |
286 |
> |
int kk = 0; |
287 |
> |
int ll = 0; |
288 |
> |
int ILscore; |
289 |
> |
int JKscore; |
290 |
> |
|
291 |
> |
std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions; |
292 |
> |
|
293 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
294 |
> |
kk = 0; |
295 |
> |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
296 |
> |
ii = 0; |
297 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
298 |
> |
ll = 0; |
299 |
> |
for (l = at4Chain.begin(); l != at4Chain.end(); l++) { |
300 |
> |
|
301 |
> |
ILscore = ii + ll; |
302 |
> |
JKscore = jj + kk; |
303 |
> |
|
304 |
> |
std::vector<std::string> myKeys; |
305 |
> |
myKeys.push_back((*i)->getName()); |
306 |
> |
myKeys.push_back((*j)->getName()); |
307 |
> |
myKeys.push_back((*k)->getName()); |
308 |
> |
myKeys.push_back((*l)->getName()); |
309 |
> |
|
310 |
> |
TorsionType* torsionType = torsionTypeCont_.find(myKeys); |
311 |
> |
if (torsionType) { |
312 |
> |
foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) ); |
313 |
> |
} |
314 |
> |
ll++; |
315 |
> |
} |
316 |
> |
ii++; |
317 |
> |
} |
318 |
> |
kk++; |
319 |
> |
} |
320 |
> |
jj++; |
321 |
> |
} |
322 |
> |
|
323 |
> |
if (foundTorsions.size() > 0) { |
324 |
> |
std::sort(foundTorsions.begin(), foundTorsions.end()); |
325 |
> |
int jkscore = foundTorsions[0].first; |
326 |
> |
int ilscore = foundTorsions[0].second; |
327 |
> |
std::vector<std::string> theKeys = foundTorsions[0].third; |
328 |
> |
|
329 |
> |
TorsionType* bestType = torsionTypeCont_.find(theKeys); |
330 |
> |
return bestType; |
331 |
> |
} else { |
332 |
> |
//if no exact match found, try wild card match |
333 |
> |
return torsionTypeCont_.find(keys, wildCardAtomTypeName_); |
334 |
> |
} |
335 |
> |
} |
336 |
> |
} |
337 |
> |
|
338 |
> |
InversionType* ForceField::getInversionType(const std::string &at1, |
339 |
> |
const std::string &at2, |
340 |
> |
const std::string &at3, |
341 |
> |
const std::string &at4) { |
342 |
> |
std::vector<std::string> keys; |
343 |
> |
keys.push_back(at1); |
344 |
> |
keys.push_back(at2); |
345 |
> |
keys.push_back(at3); |
346 |
> |
keys.push_back(at4); |
347 |
> |
|
348 |
> |
//try exact match first |
349 |
> |
InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys); |
350 |
> |
if (inversionType) { |
351 |
> |
return inversionType; |
352 |
> |
} else { |
353 |
> |
|
354 |
> |
AtomType* atype1; |
355 |
> |
AtomType* atype2; |
356 |
> |
AtomType* atype3; |
357 |
> |
AtomType* atype4; |
358 |
> |
std::vector<std::string> at1key; |
359 |
> |
at1key.push_back(at1); |
360 |
> |
atype1 = atomTypeCont_.find(at1key); |
361 |
> |
|
362 |
> |
std::vector<std::string> at2key; |
363 |
> |
at2key.push_back(at2); |
364 |
> |
atype2 = atomTypeCont_.find(at2key); |
365 |
> |
|
366 |
> |
std::vector<std::string> at3key; |
367 |
> |
at3key.push_back(at3); |
368 |
> |
atype3 = atomTypeCont_.find(at3key); |
369 |
> |
|
370 |
> |
std::vector<std::string> at4key; |
371 |
> |
at4key.push_back(at4); |
372 |
> |
atype4 = atomTypeCont_.find(at4key); |
373 |
> |
|
374 |
> |
// query atom types for their chains of responsibility |
375 |
> |
std::vector<AtomType*> at1Chain = atype1->allYourBase(); |
376 |
> |
std::vector<AtomType*> at2Chain = atype2->allYourBase(); |
377 |
> |
std::vector<AtomType*> at3Chain = atype3->allYourBase(); |
378 |
> |
std::vector<AtomType*> at4Chain = atype4->allYourBase(); |
379 |
> |
|
380 |
> |
std::vector<AtomType*>::iterator i; |
381 |
> |
std::vector<AtomType*>::iterator j; |
382 |
> |
std::vector<AtomType*>::iterator k; |
383 |
> |
std::vector<AtomType*>::iterator l; |
384 |
> |
|
385 |
> |
int ii = 0; |
386 |
> |
int jj = 0; |
387 |
> |
int kk = 0; |
388 |
> |
int ll = 0; |
389 |
> |
int Iscore; |
390 |
> |
int JKLscore; |
391 |
> |
|
392 |
> |
std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions; |
393 |
> |
|
394 |
> |
for (j = at2Chain.begin(); j != at2Chain.end(); j++) { |
395 |
> |
kk = 0; |
396 |
> |
for (k = at3Chain.begin(); k != at3Chain.end(); k++) { |
397 |
> |
ii = 0; |
398 |
> |
for (i = at1Chain.begin(); i != at1Chain.end(); i++) { |
399 |
> |
ll = 0; |
400 |
> |
for (l = at4Chain.begin(); l != at4Chain.end(); l++) { |
401 |
> |
|
402 |
> |
Iscore = ii; |
403 |
> |
JKLscore = jj + kk + ll; |
404 |
> |
|
405 |
> |
std::vector<std::string> myKeys; |
406 |
> |
myKeys.push_back((*i)->getName()); |
407 |
> |
myKeys.push_back((*j)->getName()); |
408 |
> |
myKeys.push_back((*k)->getName()); |
409 |
> |
myKeys.push_back((*l)->getName()); |
410 |
> |
|
411 |
> |
InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys); |
412 |
> |
if (inversionType) { |
413 |
> |
foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) ); |
414 |
> |
} |
415 |
> |
ll++; |
416 |
> |
} |
417 |
> |
ii++; |
418 |
> |
} |
419 |
> |
kk++; |
420 |
> |
} |
421 |
> |
jj++; |
422 |
> |
} |
423 |
> |
|
424 |
> |
if (foundInversions.size() > 0) { |
425 |
> |
std::sort(foundInversions.begin(), foundInversions.end()); |
426 |
> |
int iscore = foundInversions[0].first; |
427 |
> |
int jklscore = foundInversions[0].second; |
428 |
> |
std::vector<std::string> theKeys = foundInversions[0].third; |
429 |
> |
|
430 |
> |
InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys); |
431 |
> |
return bestType; |
432 |
> |
} else { |
433 |
> |
//if no exact match found, try wild card match |
434 |
> |
return inversionTypeCont_.find(keys, wildCardAtomTypeName_); |
435 |
> |
} |
436 |
> |
} |
437 |
> |
} |
438 |
> |
|
439 |
> |
NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) { |
440 |
> |
std::vector<std::string> keys; |
441 |
> |
keys.push_back(at1); |
442 |
> |
keys.push_back(at2); |
443 |
> |
|
444 |
> |
//try exact match first |
445 |
> |
NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys); |
446 |
> |
if (nbiType) { |
447 |
> |
return nbiType; |
448 |
> |
} else { |
449 |
> |
//if no exact match found, try wild card match |
450 |
> |
return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_); |
451 |
> |
} |
452 |
> |
} |
453 |
> |
|
454 |
> |
BondType* ForceField::getExactBondType(const std::string &at1, |
455 |
> |
const std::string &at2){ |
456 |
> |
std::vector<std::string> keys; |
457 |
> |
keys.push_back(at1); |
458 |
> |
keys.push_back(at2); |
459 |
> |
return bondTypeCont_.find(keys); |
460 |
> |
} |
461 |
> |
|
462 |
> |
BendType* ForceField::getExactBendType(const std::string &at1, |
463 |
> |
const std::string &at2, |
464 |
> |
const std::string &at3){ |
465 |
> |
std::vector<std::string> keys; |
466 |
> |
keys.push_back(at1); |
467 |
> |
keys.push_back(at2); |
468 |
> |
keys.push_back(at3); |
469 |
> |
return bendTypeCont_.find(keys); |
470 |
> |
} |
471 |
> |
|
472 |
> |
TorsionType* ForceField::getExactTorsionType(const std::string &at1, |
473 |
> |
const std::string &at2, |
474 |
> |
const std::string &at3, |
475 |
> |
const std::string &at4){ |
476 |
> |
std::vector<std::string> keys; |
477 |
> |
keys.push_back(at1); |
478 |
> |
keys.push_back(at2); |
479 |
> |
keys.push_back(at3); |
480 |
> |
keys.push_back(at4); |
481 |
> |
return torsionTypeCont_.find(keys); |
482 |
> |
} |
483 |
> |
|
484 |
> |
InversionType* ForceField::getExactInversionType(const std::string &at1, |
485 |
> |
const std::string &at2, |
486 |
> |
const std::string &at3, |
487 |
> |
const std::string &at4){ |
488 |
> |
std::vector<std::string> keys; |
489 |
> |
keys.push_back(at1); |
490 |
> |
keys.push_back(at2); |
491 |
> |
keys.push_back(at3); |
492 |
> |
keys.push_back(at4); |
493 |
> |
return inversionTypeCont_.find(keys); |
494 |
> |
} |
495 |
> |
|
496 |
> |
NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){ |
497 |
> |
std::vector<std::string> keys; |
498 |
> |
keys.push_back(at1); |
499 |
> |
keys.push_back(at2); |
500 |
> |
return nonBondedInteractionTypeCont_.find(keys); |
501 |
> |
} |
502 |
> |
|
503 |
> |
|
504 |
> |
bool ForceField::addAtomType(const std::string &at, AtomType* atomType) { |
505 |
> |
std::vector<std::string> keys; |
506 |
> |
keys.push_back(at); |
507 |
> |
return atomTypeCont_.add(keys, atomType); |
508 |
> |
} |
509 |
> |
|
510 |
> |
bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) { |
511 |
> |
std::vector<std::string> keys; |
512 |
> |
keys.push_back(at); |
513 |
> |
return atomTypeCont_.replace(keys, atomType); |
514 |
> |
} |
515 |
> |
|
516 |
> |
bool ForceField::addBondType(const std::string &at1, const std::string &at2, |
517 |
> |
BondType* bondType) { |
518 |
> |
std::vector<std::string> keys; |
519 |
> |
keys.push_back(at1); |
520 |
> |
keys.push_back(at2); |
521 |
> |
return bondTypeCont_.add(keys, bondType); |
522 |
> |
} |
523 |
> |
|
524 |
> |
bool ForceField::addBendType(const std::string &at1, const std::string &at2, |
525 |
> |
const std::string &at3, BendType* bendType) { |
526 |
> |
std::vector<std::string> keys; |
527 |
> |
keys.push_back(at1); |
528 |
> |
keys.push_back(at2); |
529 |
> |
keys.push_back(at3); |
530 |
> |
return bendTypeCont_.add(keys, bendType); |
531 |
> |
} |
532 |
> |
|
533 |
> |
bool ForceField::addTorsionType(const std::string &at1, |
534 |
> |
const std::string &at2, |
535 |
> |
const std::string &at3, |
536 |
> |
const std::string &at4, |
537 |
> |
TorsionType* torsionType) { |
538 |
> |
std::vector<std::string> keys; |
539 |
> |
keys.push_back(at1); |
540 |
> |
keys.push_back(at2); |
541 |
> |
keys.push_back(at3); |
542 |
> |
keys.push_back(at4); |
543 |
> |
return torsionTypeCont_.add(keys, torsionType); |
544 |
> |
} |
545 |
> |
|
546 |
> |
bool ForceField::addInversionType(const std::string &at1, |
547 |
> |
const std::string &at2, |
548 |
> |
const std::string &at3, |
549 |
> |
const std::string &at4, |
550 |
> |
InversionType* inversionType) { |
551 |
> |
std::vector<std::string> keys; |
552 |
> |
keys.push_back(at1); |
553 |
> |
keys.push_back(at2); |
554 |
> |
keys.push_back(at3); |
555 |
> |
keys.push_back(at4); |
556 |
> |
return inversionTypeCont_.add(keys, inversionType); |
557 |
> |
} |
558 |
> |
|
559 |
> |
bool ForceField::addNonBondedInteractionType(const std::string &at1, |
560 |
> |
const std::string &at2, |
561 |
> |
NonBondedInteractionType* nbiType) { |
562 |
> |
std::vector<std::string> keys; |
563 |
> |
keys.push_back(at1); |
564 |
> |
keys.push_back(at2); |
565 |
> |
return nonBondedInteractionTypeCont_.add(keys, nbiType); |
566 |
> |
} |
567 |
> |
|
568 |
> |
RealType ForceField::getRcutFromAtomType(AtomType* at) { |
569 |
> |
/**@todo */ |
570 |
> |
GenericData* data; |
571 |
> |
RealType rcut = 0.0; |
572 |
> |
|
573 |
> |
if (at->isLennardJones()) { |
574 |
> |
data = at->getPropertyByName("LennardJones"); |
575 |
> |
if (data != NULL) { |
576 |
> |
LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data); |
577 |
> |
|
578 |
> |
if (ljData != NULL) { |
579 |
> |
LJParam ljParam = ljData->getData(); |
580 |
> |
|
581 |
> |
//by default use 2.5*sigma as cutoff radius |
582 |
> |
rcut = 2.5 * ljParam.sigma; |
583 |
> |
|
584 |
> |
} else { |
585 |
> |
sprintf( painCave.errMsg, |
586 |
> |
"Can not cast GenericData to LJParam\n"); |
587 |
> |
painCave.severity = OPENMD_ERROR; |
588 |
> |
painCave.isFatal = 1; |
589 |
> |
simError(); |
590 |
> |
} |
591 |
> |
} else { |
592 |
> |
sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n"); |
593 |
> |
painCave.severity = OPENMD_ERROR; |
594 |
> |
painCave.isFatal = 1; |
595 |
> |
simError(); |
596 |
> |
} |
597 |
> |
} |
598 |
> |
return rcut; |
599 |
> |
} |
600 |
> |
|
601 |
> |
|
602 |
> |
ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { |
603 |
> |
std::string forceFieldFilename(filename); |
604 |
> |
ifstrstream* ffStream = new ifstrstream(); |
605 |
> |
|
606 |
> |
//try to open the force filed file in current directory first |
607 |
> |
ffStream->open(forceFieldFilename.c_str()); |
608 |
> |
if(!ffStream->is_open()){ |
609 |
> |
|
610 |
> |
forceFieldFilename = ffPath_ + "/" + forceFieldFilename; |
611 |
> |
ffStream->open( forceFieldFilename.c_str() ); |
612 |
> |
|
613 |
> |
//if current directory does not contain the force field file, |
614 |
> |
//try to open it in the path |
615 |
> |
if(!ffStream->is_open()){ |
616 |
> |
|
617 |
> |
sprintf( painCave.errMsg, |
618 |
> |
"Error opening the force field parameter file:\n" |
619 |
> |
"\t%s\n" |
620 |
> |
"\tHave you tried setting the FORCE_PARAM_PATH environment " |
621 |
> |
"variable?\n", |
622 |
> |
forceFieldFilename.c_str() ); |
623 |
> |
painCave.severity = OPENMD_ERROR; |
624 |
> |
painCave.isFatal = 1; |
625 |
> |
simError(); |
626 |
> |
} |
627 |
> |
} |
628 |
> |
return ffStream; |
629 |
> |
} |
630 |
> |
|
631 |
> |
} //end namespace OpenMD |