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] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). |
40 |
+ |
* [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). |
41 |
+ |
*/ |
42 |
+ |
|
43 |
+ |
/** |
44 |
+ |
* @file ForceField.cpp |
45 |
+ |
* @author tlin |
46 |
+ |
* @date 11/04/2004 |
47 |
+ |
* @time 22:51am |
48 |
+ |
* @version 1.0 |
49 |
+ |
*/ |
50 |
+ |
|
51 |
+ |
#include <algorithm> |
52 |
|
#include "UseTheForce/ForceField.hpp" |
53 |
+ |
#include "utils/simError.h" |
54 |
+ |
#include "utils/Tuple.hpp" |
55 |
+ |
#include "types/LennardJonesAdapter.hpp" |
56 |
|
|
57 |
< |
AtomType* ForceField::getMatchingAtomType(const string &at) { |
57 |
> |
namespace OpenMD { |
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 |
> |
ForceField::ForceField() { |
60 |
> |
|
61 |
> |
char* tempPath; |
62 |
> |
tempPath = getenv("FORCE_PARAM_PATH"); |
63 |
> |
|
64 |
> |
if (tempPath == NULL) { |
65 |
> |
//convert a macro from compiler to a string in c++ |
66 |
> |
STR_DEFINE(ffPath_, FRC_PATH ); |
67 |
> |
} else { |
68 |
> |
ffPath_ = tempPath; |
69 |
> |
} |
70 |
|
} |
13 |
– |
} |
71 |
|
|
72 |
< |
BondType* ForceField::getMatchingBondType(const string &at1, |
73 |
< |
const string &at2) { |
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 |
< |
map<pair<string,string>, BondType*>::iterator iter; |
85 |
< |
vector<BondType*> foundTypes; |
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 |
< |
iter = bondTypeMap.find(pair<at1, at2>); |
96 |
< |
if (iter != bondTypeMap.end()) { |
97 |
< |
// exact match, so just return it |
98 |
< |
return iter->second; |
99 |
< |
} |
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); |
100 |
|
|
101 |
< |
iter = bondTypeMap.find(pair<at2, at1>); |
102 |
< |
if (iter != bondTypeMap.end()) { |
103 |
< |
// exact match in reverse order, so just return it |
104 |
< |
return iter->second; |
105 |
< |
} |
101 |
> |
//try exact match first |
102 |
> |
BondType* bondType = bondTypeCont_.find(keys); |
103 |
> |
if (bondType) { |
104 |
> |
return bondType; |
105 |
> |
} else { |
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 |
< |
iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>); |
117 |
< |
if (iter != bondTypeMap.end()) { |
118 |
< |
foundTypes.push_back(iter->second); |
36 |
< |
} |
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 |
< |
iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>); |
121 |
< |
if (iter != bondTypeMap.end()) { |
40 |
< |
foundTypes.push_back(iter->second); |
41 |
< |
} |
120 |
> |
std::vector<AtomType*>::iterator i; |
121 |
> |
std::vector<AtomType*>::iterator j; |
122 |
|
|
123 |
< |
iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>); |
124 |
< |
if (iter != bondTypeMap.end()) { |
125 |
< |
foundTypes.push_back(iter->second); |
46 |
< |
} |
123 |
> |
int ii = 0; |
124 |
> |
int jj = 0; |
125 |
> |
int bondTypeScore; |
126 |
|
|
127 |
< |
iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>); |
128 |
< |
if (iter != bondTypeMap.end()) { |
129 |
< |
foundTypes.push_back(iter->second); |
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 |
< |
if (foundTypes.empty()) { |
167 |
< |
return NULL; |
168 |
< |
} else { |
169 |
< |
|
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); |
171 |
> |
keys.push_back(at2); |
172 |
> |
keys.push_back(at3); |
173 |
|
|
174 |
< |
|
174 |
> |
//try exact match first |
175 |
> |
BendType* bendType = bendTypeCont_.find(keys); |
176 |
> |
if (bendType) { |
177 |
> |
return bendType; |
178 |
> |
} else { |
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, |
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 { |
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 |
< |
BendType* ForceField::getMatchingBendType(const string &at1, const string &at2, |
284 |
< |
const string &at3); |
285 |
< |
TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2, |
68 |
< |
const string &at3, const string &at4); |
283 |
> |
std::vector<std::string> at4key; |
284 |
> |
at4key.push_back(at4); |
285 |
> |
atype4 = atomTypeCont_.find(at4key); |
286 |
|
|
287 |
< |
double ForceField::getRcutForAtomType(AtomType* at); |
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 |
< |
vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) { |
299 |
< |
|
300 |
< |
vector<vector<string> > results; |
298 |
> |
int ii = 0; |
299 |
> |
int jj = 0; |
300 |
> |
int kk = 0; |
301 |
> |
int ll = 0; |
302 |
> |
int ILscore; |
303 |
> |
int JKscore; |
304 |
|
|
305 |
< |
|
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 |
< |
vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) { |
319 |
< |
|
320 |
< |
int nStrings; |
321 |
< |
vector<string> oneResult; |
322 |
< |
vector<vector<string> > allResults; |
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 |
< |
nStrings = myAts.size(); |
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 |
< |
if (nStrings == 1) { |
353 |
< |
oneResult.push_back(wildcardCharacter); |
354 |
< |
allResults.push_back(oneResult); |
355 |
< |
return allResults; |
356 |
< |
} else { |
357 |
< |
|
358 |
< |
for (i=0; i < nStrings; i++) { |
359 |
< |
oneResult = myAts; |
360 |
< |
replace(oneResult.begin(), oneResult.end(), |
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, |
532 |
> |
const std::string &at2, |
533 |
> |
const std::string &at3){ |
534 |
> |
std::vector<std::string> keys; |
535 |
> |
keys.push_back(at1); |
536 |
> |
keys.push_back(at2); |
537 |
> |
keys.push_back(at3); |
538 |
> |
return bendTypeCont_.find(keys); |
539 |
> |
} |
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); |
548 |
> |
keys.push_back(at3); |
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::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); |
593 |
> |
} |
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; |
598 |
> |
keys.push_back(at1); |
599 |
> |
keys.push_back(at2); |
600 |
> |
keys.push_back(at3); |
601 |
> |
return bendTypeCont_.add(keys, bendType); |
602 |
> |
} |
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); |
612 |
> |
keys.push_back(at3); |
613 |
> |
keys.push_back(at4); |
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) { |
640 |
> |
RealType rcut = 0.0; |
641 |
> |
|
642 |
> |
LennardJonesAdapter lja = LennardJonesAdapter(at); |
643 |
> |
if (lja.isLennardJones()) { |
644 |
> |
rcut = 2.5 * lja.getSigma(); |
645 |
> |
} |
646 |
> |
return rcut; |
647 |
> |
} |
648 |
> |
|
649 |
> |
|
650 |
> |
ifstrstream* ForceField::openForceFieldFile(const std::string& filename) { |
651 |
> |
std::string forceFieldFilename(filename); |
652 |
> |
ifstrstream* ffStream = new ifstrstream(); |
653 |
> |
|
654 |
> |
//try to open the force filed file in current directory first |
655 |
> |
ffStream->open(forceFieldFilename.c_str()); |
656 |
> |
if(!ffStream->is_open()){ |
657 |
> |
|
658 |
> |
forceFieldFilename = ffPath_ + "/" + forceFieldFilename; |
659 |
> |
ffStream->open( forceFieldFilename.c_str() ); |
660 |
> |
|
661 |
> |
//if current directory does not contain the force field file, |
662 |
> |
//try to open it in the path |
663 |
> |
if(!ffStream->is_open()){ |
664 |
> |
|
665 |
> |
sprintf( painCave.errMsg, |
666 |
> |
"Error opening the force field parameter file:\n" |
667 |
> |
"\t%s\n" |
668 |
> |
"\tHave you tried setting the FORCE_PARAM_PATH environment " |
669 |
> |
"variable?\n", |
670 |
> |
forceFieldFilename.c_str() ); |
671 |
> |
painCave.severity = OPENMD_ERROR; |
672 |
> |
painCave.isFatal = 1; |
673 |
> |
simError(); |
674 |
> |
} |
675 |
> |
} |
676 |
> |
return ffStream; |
677 |
> |
} |
678 |
> |
|
679 |
> |
} //end namespace OpenMD |