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 1269 by gezelter, Tue Jul 1 13:28:23 2008 UTC vs.
Revision 1390 by gezelter, Wed Nov 25 20:02:06 2009 UTC

# Line 6 | Line 6
6   * redistribute this software in source and binary code form, provided
7   * that the following conditions are met:
8   *
9 < * 1. Acknowledgement of the program authors must be made in any
10 < *    publication of scientific results based in part on use of the
11 < *    program.  An acceptable form of acknowledgement is citation of
12 < *    the article in which the program was described (Matthew
13 < *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 < *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 < *    Parallel Simulation Engine for Molecular Dynamics,"
16 < *    J. Comput. Chem. 26, pp. 252-271 (2005))
17 < *
18 < * 2. Redistributions of source code must retain the above copyright
9 > * 1. Redistributions of source code must retain the above copyright
10   *    notice, this list of conditions and the following disclaimer.
11   *
12 < * 3. Redistributions in binary form must reproduce the above copyright
12 > * 2. Redistributions in binary form must reproduce the above copyright
13   *    notice, this list of conditions and the following disclaimer in the
14   *    documentation and/or other materials provided with the
15   *    distribution.
# Line 37 | Line 28
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   /**
# Line 54 | Line 54
54   #include "UseTheForce/DarkSide/atype_interface.h"
55   #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56   #include "UseTheForce/DarkSide/switcheroo_interface.h"
57 < namespace oopse {
57 > namespace OpenMD {
58  
59    ForceField::ForceField() {
60      char* tempPath;
# Line 133 | Line 133 | namespace oopse {
133          ii++;
134        }
135  
136 –      // sort the foundBonds by the score:
136  
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 <      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
143 <      BondType* bestType = bondTypeCont_.find(theKeys);
144 <      if (bestType)
145 <        return bestType;
146 <      else {
147 <        //if no exact match found, try wild card match
148 <        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
137 >      if (foundBonds.size() > 0) {
138 >        // sort the foundBonds by the score:
139 >        std::sort(foundBonds.begin(), foundBonds.end());
140 >    
141 >        int bestScore = foundBonds[0].first;
142 >        std::vector<std::string> theKeys = foundBonds[0].second;
143 >        
144 >        BondType* bestType = bondTypeCont_.find(theKeys);
145 >        
146 >        return bestType;
147 >      } else {
148 >        //if no exact match found, try wild card match
149 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
150        }
151      }
152    }
# Line 220 | Line 220 | namespace oopse {
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 <      
223 >      if (foundBends.size() > 0) {
224 >        std::sort(foundBends.begin(), foundBends.end());
225 >        int jscore = foundBends[0].first;
226 >        int ikscore = foundBends[0].second;
227 >        std::vector<std::string> theKeys = foundBends[0].third;      
228 >        
229 >        BendType* bestType = bendTypeCont_.find(theKeys);  
230 >        return bestType;
231 >      } else {        
232          //if no exact match found, try wild card match
233          return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
234        }
235      }
236    }
237  
242 –
238    TorsionType* ForceField::getTorsionType(const std::string &at1,
239                                            const std::string &at2,
240                                            const std::string &at3,
# Line 327 | Line 322 | namespace oopse {
322          jj++;
323        }
324        
325 <      std::sort(foundTorsions.begin(), foundTorsions.end());
325 >      if (foundTorsions.size() > 0) {
326 >        std::sort(foundTorsions.begin(), foundTorsions.end());
327 >        int jkscore = foundTorsions[0].first;
328 >        int ilscore = foundTorsions[0].second;
329 >        std::vector<std::string> theKeys = foundTorsions[0].third;
330 >        
331 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
332 >        return bestType;
333 >      } else {
334 >        //if no exact match found, try wild card match
335 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
336 >      }
337 >    }
338 >  }
339  
340 <      int jkscore = foundTorsions[0].first;
341 <      int ilscore = foundTorsions[0].second;
342 <      std::vector<std::string> theKeys = foundTorsions[0].third;
340 >  InversionType* ForceField::getInversionType(const std::string &at1,
341 >                                              const std::string &at2,
342 >                                              const std::string &at3,
343 >                                              const std::string &at4) {
344 >    std::vector<std::string> keys;
345 >    keys.push_back(at1);
346 >    keys.push_back(at2);    
347 >    keys.push_back(at3);    
348 >    keys.push_back(at4);    
349  
350 <      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
350 >    //try exact match first
351 >    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
352 >    if (inversionType) {
353 >      return inversionType;
354 >    } else {
355 >      
356 >      AtomType* atype1;
357 >      AtomType* atype2;
358 >      AtomType* atype3;
359 >      AtomType* atype4;
360 >      std::vector<std::string> at1key;
361 >      at1key.push_back(at1);
362 >      atype1 = atomTypeCont_.find(at1key);
363 >      
364 >      std::vector<std::string> at2key;
365 >      at2key.push_back(at2);
366 >      atype2 = atomTypeCont_.find(at2key);
367 >      
368 >      std::vector<std::string> at3key;
369 >      at3key.push_back(at3);
370 >      atype3 = atomTypeCont_.find(at3key);
371 >      
372 >      std::vector<std::string> at4key;
373 >      at4key.push_back(at4);
374 >      atype4 = atomTypeCont_.find(at4key);
375  
376 +      // query atom types for their chains of responsibility
377 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
378 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
379 +      std::vector<AtomType*> at3Chain = atype3->allYourBase();
380 +      std::vector<AtomType*> at4Chain = atype4->allYourBase();
381 +
382 +      std::vector<AtomType*>::iterator i;
383 +      std::vector<AtomType*>::iterator j;
384 +      std::vector<AtomType*>::iterator k;
385 +      std::vector<AtomType*>::iterator l;
386 +
387 +      int ii = 0;
388 +      int jj = 0;
389 +      int kk = 0;
390 +      int ll = 0;
391 +      int Iscore;
392 +      int JKLscore;
393        
394 <      TorsionType* bestType = torsionTypeCont_.find(theKeys);
395 <      if (bestType) {
396 <        return bestType;
394 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
395 >      
396 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
397 >        kk = 0;
398 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
399 >          ii = 0;      
400 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
401 >            ll = 0;
402 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
403 >              
404 >              Iscore = ii;
405 >              JKLscore = jj + kk + ll;
406 >              
407 >              std::vector<std::string> myKeys;
408 >              myKeys.push_back((*i)->getName());
409 >              myKeys.push_back((*j)->getName());
410 >              myKeys.push_back((*k)->getName());
411 >              myKeys.push_back((*l)->getName());
412 >              
413 >              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
414 >              if (inversionType) {
415 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
416 >              }
417 >              ll++;
418 >            }
419 >            ii++;
420 >          }
421 >          kk++;
422 >        }
423 >        jj++;
424 >      }
425 >        
426 >      if (foundInversions.size() > 0) {
427 >        std::sort(foundInversions.begin(), foundInversions.end());
428 >        int iscore = foundInversions[0].first;
429 >        int jklscore = foundInversions[0].second;
430 >        std::vector<std::string> theKeys = foundInversions[0].third;
431 >        
432 >        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
433 >        return bestType;
434        } else {
435          //if no exact match found, try wild card match
436 <        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
436 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
437        }
438      }
439    }
440 <
440 >  
441    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
442      std::vector<std::string> keys;
443      keys.push_back(at1);
# Line 390 | Line 482 | namespace oopse {
482      keys.push_back(at4);  
483      return torsionTypeCont_.find(keys);
484    }
485 <
485 >  
486 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
487 >                                                   const std::string &at2,
488 >                                                   const std::string &at3,
489 >                                                   const std::string &at4){
490 >    std::vector<std::string> keys;
491 >    keys.push_back(at1);
492 >    keys.push_back(at2);    
493 >    keys.push_back(at3);    
494 >    keys.push_back(at4);  
495 >    return inversionTypeCont_.find(keys);
496 >  }
497 >  
498    NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
499      std::vector<std::string> keys;
500      keys.push_back(at1);
501      keys.push_back(at2);    
502      return nonBondedInteractionTypeCont_.find(keys);
503    }
504 +  
505  
401 –
506    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
507      std::vector<std::string> keys;
508      keys.push_back(at);
509      return atomTypeCont_.add(keys, atomType);
510    }
511  
512 +  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
513 +    std::vector<std::string> keys;
514 +    keys.push_back(at);
515 +    return atomTypeCont_.replace(keys, atomType);
516 +  }
517 +
518    bool ForceField::addBondType(const std::string &at1, const std::string &at2,
519                                 BondType* bondType) {
520      std::vector<std::string> keys;
# Line 435 | Line 545 | namespace oopse {
545      return torsionTypeCont_.add(keys, torsionType);
546    }
547  
548 +  bool ForceField::addInversionType(const std::string &at1,
549 +                                    const std::string &at2,
550 +                                    const std::string &at3,
551 +                                    const std::string &at4,
552 +                                    InversionType* inversionType) {
553 +    std::vector<std::string> keys;
554 +    keys.push_back(at1);
555 +    keys.push_back(at2);    
556 +    keys.push_back(at3);    
557 +    keys.push_back(at4);    
558 +    return inversionTypeCont_.add(keys, inversionType);
559 +  }
560 +  
561    bool ForceField::addNonBondedInteractionType(const std::string &at1,
562                                                 const std::string &at2,
563                                                 NonBondedInteractionType* nbiType) {
# Line 463 | Line 586 | namespace oopse {
586          } else {
587            sprintf( painCave.errMsg,
588                     "Can not cast GenericData to LJParam\n");
589 <          painCave.severity = OOPSE_ERROR;
589 >          painCave.severity = OPENMD_ERROR;
590            painCave.isFatal = 1;
591            simError();          
592          }            
593        } else {
594          sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
595 <        painCave.severity = OOPSE_ERROR;
595 >        painCave.severity = OPENMD_ERROR;
596          painCave.isFatal = 1;
597          simError();          
598        }
# Line 499 | Line 622 | namespace oopse {
622                   "\tHave you tried setting the FORCE_PARAM_PATH environment "
623                   "variable?\n",
624                   forceFieldFilename.c_str() );
625 <        painCave.severity = OOPSE_ERROR;
625 >        painCave.severity = OPENMD_ERROR;
626          painCave.isFatal = 1;
627          simError();
628        }
# Line 512 | Line 635 | namespace oopse {
635      forceFieldOptions_.makeFortranOptions(theseFortranOptions);
636      setfForceOptions(&theseFortranOptions);
637    }
638 < } //end namespace oopse
638 > } //end namespace OpenMD

Diff Legend

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