ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/brains/ForceField.cpp
(Generate patch)

Comparing trunk/src/UseTheForce/ForceField.cpp (file contents):
Revision 1195 by cpuglis, Thu Dec 6 20:04:02 2007 UTC vs.
Revision 1277 by gezelter, Mon Jul 14 12:35:58 2008 UTC

# Line 47 | Line 47
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   #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56   #include "UseTheForce/DarkSide/switcheroo_interface.h"
# Line 89 | Line 91 | namespace oopse {
91      if (bondType) {
92        return bondType;
93      } else {
94 <      //if no exact match found, try wild card match
95 <      return bondTypeCont_.find(keys, wildCardAtomTypeName_);
94 >      AtomType* atype1;
95 >      AtomType* atype2;
96 >      std::vector<std::string> at1key;
97 >      at1key.push_back(at1);
98 >      atype1 = atomTypeCont_.find(at1key);
99 >  
100 >      std::vector<std::string> at2key;
101 >      at2key.push_back(at2);
102 >      atype2 = atomTypeCont_.find(at2key);
103 >
104 >      // query atom types for their chains of responsibility
105 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
106 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
107 >
108 >      std::vector<AtomType*>::iterator i;
109 >      std::vector<AtomType*>::iterator j;
110 >
111 >      int ii = 0;
112 >      int jj = 0;
113 >      int bondTypeScore;
114 >
115 >      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
116 >
117 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
118 >        jj = 0;
119 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
120 >
121 >          bondTypeScore = ii + jj;
122 >
123 >          std::vector<std::string> myKeys;
124 >          myKeys.push_back((*i)->getName());
125 >          myKeys.push_back((*j)->getName());
126 >
127 >          std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
128 >          BondType* bondType = bondTypeCont_.find(myKeys);
129 >          if (bondType) {
130 >            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
131 >          }
132 >          jj++;
133 >        }
134 >        ii++;
135 >      }
136 >
137 >
138 >      if (foundBonds.size() > 0) {
139 >        // sort the foundBonds by the score:
140 >        std::sort(foundBonds.begin(), foundBonds.end());
141 >    
142 >        int bestScore = foundBonds[0].first;
143 >        std::vector<std::string> theKeys = foundBonds[0].second;
144 >        
145 >        std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
146 >        BondType* bestType = bondTypeCont_.find(theKeys);
147 >        
148 >        return bestType;
149 >      } else {
150 >        //if no exact match found, try wild card match
151 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
152 >      }
153      }
154    }
155 <
155 >  
156    BendType* ForceField::getBendType(const std::string &at1,
157                                      const std::string &at2,
158                                      const std::string &at3) {
# Line 107 | Line 166 | namespace oopse {
166      if (bendType) {
167        return bendType;
168      } else {
169 <      //if no exact match found, try wild card match
170 <      return bendTypeCont_.find(keys, wildCardAtomTypeName_);
169 >
170 >      AtomType* atype1;
171 >      AtomType* atype2;
172 >      AtomType* atype3;
173 >      std::vector<std::string> at1key;
174 >      at1key.push_back(at1);
175 >      atype1 = atomTypeCont_.find(at1key);
176 >  
177 >      std::vector<std::string> at2key;
178 >      at2key.push_back(at2);
179 >      atype2 = atomTypeCont_.find(at2key);
180 >
181 >      std::vector<std::string> at3key;
182 >      at3key.push_back(at3);
183 >      atype3 = atomTypeCont_.find(at3key);
184 >
185 >      // query atom types for their chains of responsibility
186 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
187 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
188 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
189 >
190 >      std::vector<AtomType*>::iterator i;
191 >      std::vector<AtomType*>::iterator j;
192 >      std::vector<AtomType*>::iterator k;
193 >
194 >      int ii = 0;
195 >      int jj = 0;
196 >      int kk = 0;
197 >      int IKscore;
198 >
199 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
200 >
201 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
202 >        ii = 0;
203 >        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
204 >          kk = 0;
205 >          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
206 >          
207 >            IKscore = ii + kk;
208 >
209 >            std::vector<std::string> myKeys;
210 >            myKeys.push_back((*i)->getName());
211 >            myKeys.push_back((*j)->getName());
212 >            myKeys.push_back((*k)->getName());
213 >
214 >            BendType* bendType = bendTypeCont_.find(myKeys);
215 >            if (bendType) {
216 >              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
217 >            }
218 >            kk++;
219 >          }
220 >          ii++;
221 >        }
222 >        jj++;
223 >      }
224 >      
225 >      if (foundBends.size() > 0) {
226 >        std::sort(foundBends.begin(), foundBends.end());
227 >        int jscore = foundBends[0].first;
228 >        int ikscore = foundBends[0].second;
229 >        std::vector<std::string> theKeys = foundBends[0].third;
230 >        
231 >        std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
232 >        
233 >        BendType* bestType = bendTypeCont_.find(theKeys);  
234 >        return bestType;
235 >      } else {        
236 >        //if no exact match found, try wild card match
237 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
238 >      }
239      }
240    }
241  
# Line 122 | Line 249 | namespace oopse {
249      keys.push_back(at3);    
250      keys.push_back(at4);    
251  
252 +
253 +    //try exact match first
254      TorsionType* torsionType = torsionTypeCont_.find(keys);
255      if (torsionType) {
256        return torsionType;
257      } else {
258 <      //if no exact match found, try wild card match
259 <      return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
260 <    }
261 <    
262 <    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
258 >
259 >      AtomType* atype1;
260 >      AtomType* atype2;
261 >      AtomType* atype3;
262 >      AtomType* atype4;
263 >      std::vector<std::string> at1key;
264 >      at1key.push_back(at1);
265 >      atype1 = atomTypeCont_.find(at1key);
266 >  
267 >      std::vector<std::string> at2key;
268 >      at2key.push_back(at2);
269 >      atype2 = atomTypeCont_.find(at2key);
270 >
271 >      std::vector<std::string> at3key;
272 >      at3key.push_back(at3);
273 >      atype3 = atomTypeCont_.find(at3key);
274 >
275 >      std::vector<std::string> at4key;
276 >      at4key.push_back(at4);
277 >      atype4 = atomTypeCont_.find(at4key);
278 >
279 >      // query atom types for their chains of responsibility
280 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
281 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
282 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
283 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
284 >
285 >      std::vector<AtomType*>::iterator i;
286 >      std::vector<AtomType*>::iterator j;
287 >      std::vector<AtomType*>::iterator k;
288 >      std::vector<AtomType*>::iterator l;
289 >
290 >      int ii = 0;
291 >      int jj = 0;
292 >      int kk = 0;
293 >      int ll = 0;
294 >      int ILscore;
295 >      int JKscore;
296 >
297 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
298 >
299 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
300 >        kk = 0;
301 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
302 >          ii = 0;      
303 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
304 >            ll = 0;
305 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
306 >          
307 >              ILscore = ii + ll;
308 >              JKscore = jj + kk;
309 >
310 >              std::vector<std::string> myKeys;
311 >              myKeys.push_back((*i)->getName());
312 >              myKeys.push_back((*j)->getName());
313 >              myKeys.push_back((*k)->getName());
314 >              myKeys.push_back((*l)->getName());
315 >
316 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
317 >              if (torsionType) {
318 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
319 >              }
320 >              ll++;
321 >            }
322 >            ii++;
323 >          }
324 >          kk++;
325 >        }
326 >        jj++;
327 >      }
328 >      
329 >      if (foundTorsions.size() > 0) {
330 >        std::sort(foundTorsions.begin(), foundTorsions.end());
331 >        int jkscore = foundTorsions[0].first;
332 >        int ilscore = foundTorsions[0].second;
333 >        std::vector<std::string> theKeys = foundTorsions[0].third;
334 >        
335 >        std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
336 >                
337 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
338 >        return bestType;
339 >      } else {
340 >        //if no exact match found, try wild card match
341 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
342 >      }
343 >    }
344    }
345 +
346 +  InversionType* ForceField::getInversionType(const std::string &at1,
347 +                                              const std::string &at2,
348 +                                              const std::string &at3,
349 +                                              const std::string &at4) {
350 +    std::vector<std::string> keys;
351 +    keys.push_back(at1);
352 +    keys.push_back(at2);    
353 +    keys.push_back(at3);    
354 +    keys.push_back(at4);    
355 +
356 +    //try exact match first
357 +    InversionType* inversionType = inversionTypeCont_.find(keys);
358 +    if (inversionType) {
359 +      return inversionType;
360 +    } else {
361 +      
362 +      AtomType* atype1;
363 +      AtomType* atype2;
364 +      AtomType* atype3;
365 +      AtomType* atype4;
366 +      std::vector<std::string> at1key;
367 +      at1key.push_back(at1);
368 +      atype1 = atomTypeCont_.find(at1key);
369 +      
370 +      std::vector<std::string> at2key;
371 +      at2key.push_back(at2);
372 +      atype2 = atomTypeCont_.find(at2key);
373 +      
374 +      std::vector<std::string> at3key;
375 +      at3key.push_back(at3);
376 +      atype3 = atomTypeCont_.find(at3key);
377 +      
378 +      std::vector<std::string> at4key;
379 +      at4key.push_back(at4);
380 +      atype4 = atomTypeCont_.find(at4key);
381 +
382 +      // query atom types for their chains of responsibility
383 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
384 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
385 +      std::vector<AtomType*> at3Chain = atype3->allYourBase();
386 +      std::vector<AtomType*> at4Chain = atype4->allYourBase();
387 +
388 +      std::vector<AtomType*>::iterator i;
389 +      std::vector<AtomType*>::iterator j;
390 +      std::vector<AtomType*>::iterator k;
391 +      std::vector<AtomType*>::iterator l;
392  
393 +      int ii = 0;
394 +      int jj = 0;
395 +      int kk = 0;
396 +      int ll = 0;
397 +      int Iscore;
398 +      int JKLscore;
399 +      
400 +      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
401 +      
402 +      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
403 +        kk = 0;
404 +        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
405 +          ii = 0;      
406 +          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
407 +            ll = 0;
408 +            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
409 +              
410 +              Iscore = ii;
411 +              JKLscore = jj + kk + ll;
412 +              
413 +              std::vector<std::string> myKeys;
414 +              myKeys.push_back((*i)->getName());
415 +              myKeys.push_back((*j)->getName());
416 +              myKeys.push_back((*k)->getName());
417 +              myKeys.push_back((*l)->getName());
418 +              
419 +              InversionType* inversionType = inversionTypeCont_.find(myKeys);
420 +              if (inversionType) {
421 +                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
422 +              }
423 +              ll++;
424 +            }
425 +            ii++;
426 +          }
427 +          kk++;
428 +        }
429 +        jj++;
430 +      }
431 +        
432 +      if (foundInversions.size() > 0) {
433 +        std::sort(foundInversions.begin(), foundInversions.end());
434 +        int iscore = foundInversions[0].first;
435 +        int jklscore = foundInversions[0].second;
436 +        std::vector<std::string> theKeys = foundInversions[0].third;
437 +        
438 +        std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
439 +                
440 +        InversionType* bestType = inversionTypeCont_.find(theKeys);
441 +        return bestType;
442 +      } else {
443 +        //if no exact match found, try wild card match
444 +        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
445 +      }
446 +    }
447 +  }
448 +  
449    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
450      std::vector<std::string> keys;
451      keys.push_back(at1);
# Line 177 | Line 490 | namespace oopse {
490      keys.push_back(at4);  
491      return torsionTypeCont_.find(keys);
492    }
493 <
493 >  
494 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
495 >                                                   const std::string &at2,
496 >                                                   const std::string &at3,
497 >                                                   const std::string &at4){
498 >    std::vector<std::string> keys;
499 >    keys.push_back(at1);
500 >    keys.push_back(at2);    
501 >    keys.push_back(at3);    
502 >    keys.push_back(at4);  
503 >    return inversionTypeCont_.find(keys);
504 >  }
505 >  
506    NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
507      std::vector<std::string> keys;
508      keys.push_back(at1);
509      keys.push_back(at2);    
510      return nonBondedInteractionTypeCont_.find(keys);
511    }
512 +  
513  
188 –
514    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
515      std::vector<std::string> keys;
516      keys.push_back(at);
# Line 222 | Line 547 | namespace oopse {
547      return torsionTypeCont_.add(keys, torsionType);
548    }
549  
550 +  bool ForceField::addInversionType(const std::string &at1,
551 +                                    const std::string &at2,
552 +                                    const std::string &at3,
553 +                                    const std::string &at4,
554 +                                    InversionType* inversionType) {
555 +    std::vector<std::string> keys;
556 +    keys.push_back(at1);
557 +    keys.push_back(at2);    
558 +    keys.push_back(at3);    
559 +    keys.push_back(at4);    
560 +    return inversionTypeCont_.add(keys, inversionType);
561 +  }
562 +  
563    bool ForceField::addNonBondedInteractionType(const std::string &at1,
564                                                 const std::string &at2,
565                                                 NonBondedInteractionType* nbiType) {

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines