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 507 by gezelter, Fri Apr 15 22:04:00 2005 UTC vs.
Revision 1275 by cli2, Fri Jul 4 20:54:29 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"
57   namespace oopse {
58  
59    ForceField::ForceField() {
# Line 67 | Line 71 | namespace oopse {
71  
72    ForceField::~ForceField() {
73      deleteAtypes();
74 +    deleteSwitch();
75    }
76  
77    AtomType* ForceField::getAtomType(const std::string &at) {
# Line 75 | Line 80 | namespace oopse {
80      return atomTypeCont_.find(keys);
81    }
82  
83 <  BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
83 >  BondType* ForceField::getBondType(const std::string &at1,
84 >                                    const std::string &at2) {
85      std::vector<std::string> keys;
86      keys.push_back(at1);
87      keys.push_back(at2);    
# Line 85 | 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_);
96 <    }
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 <  }
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 <  BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
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 >          BondType* bondType = bondTypeCont_.find(myKeys);
128 >          if (bondType) {
129 >            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
130 >          }
131 >          jj++;
132 >        }
133 >        ii++;
134 >      }
135 >
136 >      // sort the foundBonds by the score:
137 >
138 >      std::sort(foundBonds.begin(), foundBonds.end());
139 >      
140 >      int bestScore = foundBonds[0].first;
141 >      std::vector<std::string> theKeys = foundBonds[0].second;
142 >
143 >      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
144 >      BondType* bestType = bondTypeCont_.find(theKeys);
145 >      if (bestType)
146 >        return bestType;
147 >      else {
148 >        //if no exact match found, try wild card match
149 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
150 >      }
151 >    }
152 >  }
153 >  
154 >  BendType* ForceField::getBendType(const std::string &at1,
155 >                                    const std::string &at2,
156                                      const std::string &at3) {
157      std::vector<std::string> keys;
158      keys.push_back(at1);
# Line 103 | Line 164 | namespace oopse {
164      if (bendType) {
165        return bendType;
166      } else {
167 <      //if no exact match found, try wild card match
168 <      return bendTypeCont_.find(keys, wildCardAtomTypeName_);
167 >
168 >      AtomType* atype1;
169 >      AtomType* atype2;
170 >      AtomType* atype3;
171 >      std::vector<std::string> at1key;
172 >      at1key.push_back(at1);
173 >      atype1 = atomTypeCont_.find(at1key);
174 >  
175 >      std::vector<std::string> at2key;
176 >      at2key.push_back(at2);
177 >      atype2 = atomTypeCont_.find(at2key);
178 >
179 >      std::vector<std::string> at3key;
180 >      at3key.push_back(at3);
181 >      atype3 = atomTypeCont_.find(at3key);
182 >
183 >      // query atom types for their chains of responsibility
184 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
185 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
186 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
187 >
188 >      std::vector<AtomType*>::iterator i;
189 >      std::vector<AtomType*>::iterator j;
190 >      std::vector<AtomType*>::iterator k;
191 >
192 >      int ii = 0;
193 >      int jj = 0;
194 >      int kk = 0;
195 >      int IKscore;
196 >
197 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
198 >
199 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
200 >        ii = 0;
201 >        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
202 >          kk = 0;
203 >          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
204 >          
205 >            IKscore = ii + kk;
206 >
207 >            std::vector<std::string> myKeys;
208 >            myKeys.push_back((*i)->getName());
209 >            myKeys.push_back((*j)->getName());
210 >            myKeys.push_back((*k)->getName());
211 >
212 >            BendType* bendType = bendTypeCont_.find(myKeys);
213 >            if (bendType) {
214 >              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
215 >            }
216 >            kk++;
217 >          }
218 >          ii++;
219 >        }
220 >        jj++;
221 >      }
222 >      
223 >      std::sort(foundBends.begin(), foundBends.end());
224 >
225 >      int jscore = foundBends[0].first;
226 >      int ikscore = foundBends[0].second;
227 >      std::vector<std::string> theKeys = foundBends[0].third;
228 >
229 >      std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
230 >
231 >      BendType* bestType = bendTypeCont_.find(theKeys);  
232 >      if (bestType)
233 >        return bestType;
234 >      else {
235 >      
236 >        //if no exact match found, try wild card match
237 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
238 >      }
239      }
240    }
241  
242 <  TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
243 <                                          const std::string &at3, const std::string &at4) {
242 >
243 >  TorsionType* ForceField::getTorsionType(const std::string &at1,
244 >                                          const std::string &at2,
245 >                                          const std::string &at3,
246 >                                          const std::string &at4) {
247      std::vector<std::string> keys;
248      keys.push_back(at1);
249      keys.push_back(at2);    
250      keys.push_back(at3);    
251      keys.push_back(at4);    
252  
253 +
254 +    //try exact match first
255      TorsionType* torsionType = torsionTypeCont_.find(keys);
256      if (torsionType) {
257        return torsionType;
258      } else {
259 <      //if no exact match found, try wild card match
260 <      return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
261 <    }
262 <    
263 <    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
259 >
260 >      AtomType* atype1;
261 >      AtomType* atype2;
262 >      AtomType* atype3;
263 >      AtomType* atype4;
264 >      std::vector<std::string> at1key;
265 >      at1key.push_back(at1);
266 >      atype1 = atomTypeCont_.find(at1key);
267 >  
268 >      std::vector<std::string> at2key;
269 >      at2key.push_back(at2);
270 >      atype2 = atomTypeCont_.find(at2key);
271 >
272 >      std::vector<std::string> at3key;
273 >      at3key.push_back(at3);
274 >      atype3 = atomTypeCont_.find(at3key);
275 >
276 >      std::vector<std::string> at4key;
277 >      at4key.push_back(at4);
278 >      atype4 = atomTypeCont_.find(at4key);
279 >
280 >      // query atom types for their chains of responsibility
281 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
282 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
283 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
284 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
285 >
286 >      std::vector<AtomType*>::iterator i;
287 >      std::vector<AtomType*>::iterator j;
288 >      std::vector<AtomType*>::iterator k;
289 >      std::vector<AtomType*>::iterator l;
290 >
291 >      int ii = 0;
292 >      int jj = 0;
293 >      int kk = 0;
294 >      int ll = 0;
295 >      int ILscore;
296 >      int JKscore;
297 >
298 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
299 >
300 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
301 >        kk = 0;
302 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
303 >          ii = 0;      
304 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
305 >            ll = 0;
306 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
307 >          
308 >              ILscore = ii + ll;
309 >              JKscore = jj + kk;
310 >
311 >              std::vector<std::string> myKeys;
312 >              myKeys.push_back((*i)->getName());
313 >              myKeys.push_back((*j)->getName());
314 >              myKeys.push_back((*k)->getName());
315 >              myKeys.push_back((*l)->getName());
316 >
317 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
318 >              if (torsionType) {
319 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
320 >              }
321 >              ll++;
322 >            }
323 >            ii++;
324 >          }
325 >          kk++;
326 >        }
327 >        jj++;
328 >      }
329 >      
330 >      std::sort(foundTorsions.begin(), foundTorsions.end());
331 >
332 >      int jkscore = foundTorsions[0].first;
333 >      int ilscore = foundTorsions[0].second;
334 >      std::vector<std::string> theKeys = foundTorsions[0].third;
335  
336 +      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
337 +
338 +      
339 +      TorsionType* bestType = torsionTypeCont_.find(theKeys);
340 +      if (bestType) {
341 +        return bestType;
342 +      } else {
343 +        //if no exact match found, try wild card match
344 +        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
345 +      }
346 +    }
347    }
348  
349 <  BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
349 >  InversionType* ForceField::getInversionType(const std::string &at1,
350 >                                              const std::string &at2,
351 >                                              const std::string &at3,
352 >                                              const std::string &at4) {
353      std::vector<std::string> keys;
354      keys.push_back(at1);
355      keys.push_back(at2);    
356 +    keys.push_back(at3);    
357 +    keys.push_back(at4);    
358 +
359 +    //try exact match first
360 +    InversionType* inversionType = inversionTypeCont_.find(keys);
361 +    if (inversionType) {
362 +      return inversionType;
363 +    } else {
364 +      
365 +      AtomType* atype1;
366 +      AtomType* atype2;
367 +      AtomType* atype3;
368 +      AtomType* atype4;
369 +      std::vector<std::string> at1key;
370 +      at1key.push_back(at1);
371 +      atype1 = atomTypeCont_.find(at1key);
372 +      
373 +      std::vector<std::string> at2key;
374 +      at2key.push_back(at2);
375 +      atype2 = atomTypeCont_.find(at2key);
376 +      
377 +      std::vector<std::string> at3key;
378 +      at3key.push_back(at3);
379 +      atype3 = atomTypeCont_.find(at3key);
380 +      
381 +      std::vector<std::string> at4key;
382 +      at4key.push_back(at4);
383 +      atype4 = atomTypeCont_.find(at4key);
384 +
385 +      // query atom types for their chains of responsibility
386 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
387 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
388 +      std::vector<AtomType*> at3Chain = atype3->allYourBase();
389 +      std::vector<AtomType*> at4Chain = atype4->allYourBase();
390 +
391 +      std::vector<AtomType*>::iterator i;
392 +      std::vector<AtomType*>::iterator j;
393 +      std::vector<AtomType*>::iterator k;
394 +      std::vector<AtomType*>::iterator l;
395 +
396 +      int ii = 0;
397 +      int jj = 0;
398 +      int kk = 0;
399 +      int ll = 0;
400 +      int Iscore;
401 +      int JKLscore;
402 +      
403 +      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
404 +      
405 +      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
406 +        kk = 0;
407 +        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
408 +          ii = 0;      
409 +          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
410 +            ll = 0;
411 +            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
412 +              
413 +              Iscore = ii;
414 +              JKLscore = jj + kk + ll;
415 +              
416 +              std::vector<std::string> myKeys;
417 +              myKeys.push_back((*i)->getName());
418 +              myKeys.push_back((*j)->getName());
419 +              myKeys.push_back((*k)->getName());
420 +              myKeys.push_back((*l)->getName());
421 +              
422 +              InversionType* inversionType = inversionTypeCont_.find(myKeys);
423 +              if (inversionType) {
424 +                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
425 +              }
426 +              ll++;
427 +            }
428 +            ii++;
429 +          }
430 +          kk++;
431 +        }
432 +        jj++;
433 +      }
434 +      
435 +      std::sort(foundInversions.begin(), foundInversions.end());
436 +      
437 +      int iscore = foundInversions[0].first;
438 +      int jklscore = foundInversions[0].second;
439 +      std::vector<std::string> theKeys = foundInversions[0].third;
440 +      
441 +      std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
442 +      
443 +      
444 +      InversionType* bestType = inversionTypeCont_.find(theKeys);
445 +      if (bestType) {
446 +        return bestType;
447 +      } else {
448 +        //if no exact match found, try wild card match
449 +        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
450 +      }
451 +    }
452 +  }
453 +  
454 +  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
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 +      //if no exact match found, try wild card match
465 +      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
466 +    }    
467 +  }
468 +  
469 +  BondType* ForceField::getExactBondType(const std::string &at1,
470 +                                         const std::string &at2){
471 +    std::vector<std::string> keys;
472 +    keys.push_back(at1);
473 +    keys.push_back(at2);    
474      return bondTypeCont_.find(keys);
475    }
476 <
477 <  BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
476 >  
477 >  BendType* ForceField::getExactBendType(const std::string &at1,
478 >                                         const std::string &at2,
479                                           const std::string &at3){
480      std::vector<std::string> keys;
481      keys.push_back(at1);
# Line 143 | Line 483 | namespace oopse {
483      keys.push_back(at3);    
484      return bendTypeCont_.find(keys);
485    }
486 <
487 <  TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
488 <                                               const std::string &at3, const std::string &at4){
486 >  
487 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
488 >                                               const std::string &at2,
489 >                                               const std::string &at3,
490 >                                               const std::string &at4){
491      std::vector<std::string> keys;
492      keys.push_back(at1);
493      keys.push_back(at2);    
# Line 153 | Line 495 | namespace oopse {
495      keys.push_back(at4);  
496      return torsionTypeCont_.find(keys);
497    }
498 +  
499 +  InversionType* ForceField::getExactInversionType(const std::string &at1,
500 +                                                   const std::string &at2,
501 +                                                   const std::string &at3,
502 +                                                   const std::string &at4){
503 +    std::vector<std::string> keys;
504 +    keys.push_back(at1);
505 +    keys.push_back(at2);    
506 +    keys.push_back(at3);    
507 +    keys.push_back(at4);  
508 +    return inversionTypeCont_.find(keys);
509 +  }
510 +  
511 +  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
512 +    std::vector<std::string> keys;
513 +    keys.push_back(at1);
514 +    keys.push_back(at2);    
515 +    return nonBondedInteractionTypeCont_.find(keys);
516 +  }
517 +  
518 +
519    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
520      std::vector<std::string> keys;
521      keys.push_back(at);
522      return atomTypeCont_.add(keys, atomType);
523    }
524  
525 <  bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
525 >  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
526 >                               BondType* bondType) {
527      std::vector<std::string> keys;
528      keys.push_back(at1);
529      keys.push_back(at2);    
530 <    return bondTypeCont_.add(keys, bondType);
167 <
530 >    return bondTypeCont_.add(keys, bondType);    
531    }
532 <
532 >  
533    bool ForceField::addBendType(const std::string &at1, const std::string &at2,
534                                 const std::string &at3, BendType* bendType) {
535      std::vector<std::string> keys;
# Line 175 | Line 538 | namespace oopse {
538      keys.push_back(at3);    
539      return bendTypeCont_.add(keys, bendType);
540    }
541 <
542 <  bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
543 <                                  const std::string &at3, const std::string &at4, TorsionType* torsionType) {
541 >  
542 >  bool ForceField::addTorsionType(const std::string &at1,
543 >                                  const std::string &at2,
544 >                                  const std::string &at3,
545 >                                  const std::string &at4,
546 >                                  TorsionType* torsionType) {
547      std::vector<std::string> keys;
548      keys.push_back(at1);
549      keys.push_back(at2);    
# Line 186 | Line 552 | namespace oopse {
552      return torsionTypeCont_.add(keys, torsionType);
553    }
554  
555 <  double ForceField::getRcutFromAtomType(AtomType* at) {
555 >  bool ForceField::addInversionType(const std::string &at1,
556 >                                    const std::string &at2,
557 >                                    const std::string &at3,
558 >                                    const std::string &at4,
559 >                                    InversionType* inversionType) {
560 >    std::vector<std::string> keys;
561 >    keys.push_back(at1);
562 >    keys.push_back(at2);    
563 >    keys.push_back(at3);    
564 >    keys.push_back(at4);    
565 >    return inversionTypeCont_.add(keys, inversionType);
566 >  }
567 >  
568 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
569 >                                               const std::string &at2,
570 >                                               NonBondedInteractionType* nbiType) {
571 >    std::vector<std::string> keys;
572 >    keys.push_back(at1);
573 >    keys.push_back(at2);    
574 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
575 >  }
576 >  
577 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
578      /**@todo */
579      GenericData* data;
580 <    double rcut = 0.0;
581 <
580 >    RealType rcut = 0.0;
581 >    
582      if (at->isLennardJones()) {
583        data = at->getPropertyByName("LennardJones");
584        if (data != NULL) {
585          LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
586 <
586 >        
587          if (ljData != NULL) {
588            LJParam ljParam = ljData->getData();
589 <
589 >          
590            //by default use 2.5*sigma as cutoff radius
591            rcut = 2.5 * ljParam.sigma;
592 <                
592 >          
593          } else {
594            sprintf( painCave.errMsg,
595                     "Can not cast GenericData to LJParam\n");
# Line 216 | Line 604 | namespace oopse {
604          simError();          
605        }
606      }
219 –
607      return rcut;    
608    }
609 +  
610  
223 –
611    ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
612      std::string forceFieldFilename(filename);
613      ifstrstream* ffStream = new ifstrstream();
# Line 247 | Line 634 | namespace oopse {
634          simError();
635        }
636      }  
250 –
637      return ffStream;
252 –
638    }
639  
640 +  void ForceField::setFortranForceOptions(){
641 +    ForceOptions theseFortranOptions;
642 +    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
643 +    setfForceOptions(&theseFortranOptions);
644 +  }
645   } //end namespace oopse

Diff Legend

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