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.
branches/development/src/UseTheForce/ForceField.cpp (file contents), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 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]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 + * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   /**
# Line 51 | Line 52
52   #include "UseTheForce/ForceField.hpp"
53   #include "utils/simError.h"
54   #include "utils/Tuple.hpp"
55 < #include "UseTheForce/DarkSide/atype_interface.h"
55 < #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56 < #include "UseTheForce/DarkSide/switcheroo_interface.h"
57 < namespace oopse {
55 > namespace OpenMD {
56  
57    ForceField::ForceField() {
58 +
59      char* tempPath;
60      tempPath = getenv("FORCE_PARAM_PATH");
61 <
61 >    
62      if (tempPath == NULL) {
63        //convert a macro from compiler to a string in c++
64        STR_DEFINE(ffPath_, FRC_PATH );
# Line 68 | Line 67 | namespace oopse {
67      }
68    }
69  
70 <
71 <  ForceField::~ForceField() {
72 <    deleteAtypes();
73 <    deleteSwitch();
74 <  }
75 <
70 >  /**
71 >   * getAtomType by string
72 >   *
73 >   * finds the requested atom type in this force field using the string
74 >   * name of the atom type.
75 >   */
76    AtomType* ForceField::getAtomType(const std::string &at) {
77      std::vector<std::string> keys;
78      keys.push_back(at);
79      return atomTypeCont_.find(keys);
80    }
81  
82 +  /**
83 +   * getAtomType by ident
84 +   *
85 +   * finds the requested atom type in this force field using the
86 +   * integer ident instead of the string name of the atom type.
87 +   */
88 +  AtomType* ForceField::getAtomType(int ident) {  
89 +    std::string at = atypeIdentToName.find(ident)->second;
90 +    return getAtomType(at);
91 +  }
92 +
93    BondType* ForceField::getBondType(const std::string &at1,
94                                      const std::string &at2) {
95      std::vector<std::string> keys;
# Line 133 | Line 143 | namespace oopse {
143          ii++;
144        }
145  
136      // sort the foundBonds by the score:
146  
147 <      std::sort(foundBonds.begin(), foundBonds.end());
148 <      
149 <      int bestScore = foundBonds[0].first;
150 <      std::vector<std::string> theKeys = foundBonds[0].second;
151 <
152 <      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
153 <      BondType* bestType = bondTypeCont_.find(theKeys);
154 <      if (bestType)
155 <        return bestType;
156 <      else {
157 <        //if no exact match found, try wild card match
158 <        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
147 >      if (foundBonds.size() > 0) {
148 >        // sort the foundBonds by the score:
149 >        std::sort(foundBonds.begin(), foundBonds.end());
150 >    
151 >        int bestScore = foundBonds[0].first;
152 >        std::vector<std::string> theKeys = foundBonds[0].second;
153 >        
154 >        BondType* bestType = bondTypeCont_.find(theKeys);
155 >        
156 >        return bestType;
157 >      } else {
158 >        //if no exact match found, try wild card match
159 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
160        }
161      }
162    }
# Line 220 | Line 230 | namespace oopse {
230          jj++;
231        }
232        
233 <      std::sort(foundBends.begin(), foundBends.end());
234 <
235 <      int jscore = foundBends[0].first;
236 <      int ikscore = foundBends[0].second;
237 <      std::vector<std::string> theKeys = foundBends[0].third;
238 <
239 <      std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
240 <
241 <      BendType* bestType = bendTypeCont_.find(theKeys);  
232 <      if (bestType)
233 <        return bestType;
234 <      else {
235 <      
233 >      if (foundBends.size() > 0) {
234 >        std::sort(foundBends.begin(), foundBends.end());
235 >        int jscore = foundBends[0].first;
236 >        int ikscore = foundBends[0].second;
237 >        std::vector<std::string> theKeys = foundBends[0].third;      
238 >        
239 >        BendType* bestType = bendTypeCont_.find(theKeys);  
240 >        return bestType;
241 >      } else {        
242          //if no exact match found, try wild card match
243          return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
244        }
245      }
246    }
247  
242
248    TorsionType* ForceField::getTorsionType(const std::string &at1,
249                                            const std::string &at2,
250                                            const std::string &at3,
# Line 327 | Line 332 | namespace oopse {
332          jj++;
333        }
334        
335 <      std::sort(foundTorsions.begin(), foundTorsions.end());
335 >      if (foundTorsions.size() > 0) {
336 >        std::sort(foundTorsions.begin(), foundTorsions.end());
337 >        int jkscore = foundTorsions[0].first;
338 >        int ilscore = foundTorsions[0].second;
339 >        std::vector<std::string> theKeys = foundTorsions[0].third;
340 >        
341 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
342 >        return bestType;
343 >      } else {
344 >        //if no exact match found, try wild card match
345 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
346 >      }
347 >    }
348 >  }
349  
350 <      int jkscore = foundTorsions[0].first;
351 <      int ilscore = foundTorsions[0].second;
352 <      std::vector<std::string> theKeys = foundTorsions[0].third;
350 >  InversionType* ForceField::getInversionType(const std::string &at1,
351 >                                              const std::string &at2,
352 >                                              const std::string &at3,
353 >                                              const std::string &at4) {
354 >    std::vector<std::string> keys;
355 >    keys.push_back(at1);
356 >    keys.push_back(at2);    
357 >    keys.push_back(at3);    
358 >    keys.push_back(at4);    
359  
360 <      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
360 >    //try exact match first
361 >    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
362 >    if (inversionType) {
363 >      return inversionType;
364 >    } else {
365 >      
366 >      AtomType* atype1;
367 >      AtomType* atype2;
368 >      AtomType* atype3;
369 >      AtomType* atype4;
370 >      std::vector<std::string> at1key;
371 >      at1key.push_back(at1);
372 >      atype1 = atomTypeCont_.find(at1key);
373 >      
374 >      std::vector<std::string> at2key;
375 >      at2key.push_back(at2);
376 >      atype2 = atomTypeCont_.find(at2key);
377 >      
378 >      std::vector<std::string> at3key;
379 >      at3key.push_back(at3);
380 >      atype3 = atomTypeCont_.find(at3key);
381 >      
382 >      std::vector<std::string> at4key;
383 >      at4key.push_back(at4);
384 >      atype4 = atomTypeCont_.find(at4key);
385  
386 +      // query atom types for their chains of responsibility
387 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
388 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
389 +      std::vector<AtomType*> at3Chain = atype3->allYourBase();
390 +      std::vector<AtomType*> at4Chain = atype4->allYourBase();
391 +
392 +      std::vector<AtomType*>::iterator i;
393 +      std::vector<AtomType*>::iterator j;
394 +      std::vector<AtomType*>::iterator k;
395 +      std::vector<AtomType*>::iterator l;
396 +
397 +      int ii = 0;
398 +      int jj = 0;
399 +      int kk = 0;
400 +      int ll = 0;
401 +      int Iscore;
402 +      int JKLscore;
403        
404 <      TorsionType* bestType = torsionTypeCont_.find(theKeys);
405 <      if (bestType) {
406 <        return bestType;
404 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
405 >      
406 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
407 >        kk = 0;
408 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
409 >          ii = 0;      
410 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
411 >            ll = 0;
412 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
413 >              
414 >              Iscore = ii;
415 >              JKLscore = jj + kk + ll;
416 >              
417 >              std::vector<std::string> myKeys;
418 >              myKeys.push_back((*i)->getName());
419 >              myKeys.push_back((*j)->getName());
420 >              myKeys.push_back((*k)->getName());
421 >              myKeys.push_back((*l)->getName());
422 >              
423 >              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
424 >              if (inversionType) {
425 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
426 >              }
427 >              ll++;
428 >            }
429 >            ii++;
430 >          }
431 >          kk++;
432 >        }
433 >        jj++;
434 >      }
435 >        
436 >      if (foundInversions.size() > 0) {
437 >        std::sort(foundInversions.begin(), foundInversions.end());
438 >        int iscore = foundInversions[0].first;
439 >        int jklscore = foundInversions[0].second;
440 >        std::vector<std::string> theKeys = foundInversions[0].third;
441 >        
442 >        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
443 >        return bestType;
444        } else {
445          //if no exact match found, try wild card match
446 <        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
446 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
447        }
448      }
449    }
450 <
450 >  
451    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
452 +    
453      std::vector<std::string> keys;
454      keys.push_back(at1);
455      keys.push_back(at2);    
# Line 356 | Line 459 | namespace oopse {
459      if (nbiType) {
460        return nbiType;
461      } else {
462 <      //if no exact match found, try wild card match
463 <      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
464 <    }    
462 >      AtomType* atype1;
463 >      AtomType* atype2;
464 >      std::vector<std::string> at1key;
465 >      at1key.push_back(at1);
466 >      atype1 = atomTypeCont_.find(at1key);
467 >      
468 >      std::vector<std::string> at2key;
469 >      at2key.push_back(at2);
470 >      atype2 = atomTypeCont_.find(at2key);
471 >      
472 >      // query atom types for their chains of responsibility
473 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
474 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
475 >      
476 >      std::vector<AtomType*>::iterator i;
477 >      std::vector<AtomType*>::iterator j;
478 >      
479 >      int ii = 0;
480 >      int jj = 0;
481 >      int nbiTypeScore;
482 >      
483 >      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
484 >      
485 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
486 >        jj = 0;
487 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
488 >          
489 >          nbiTypeScore = ii + jj;
490 >          
491 >          std::vector<std::string> myKeys;
492 >          myKeys.push_back((*i)->getName());
493 >          myKeys.push_back((*j)->getName());
494 >          
495 >          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
496 >          if (nbiType) {
497 >            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
498 >          }
499 >          jj++;
500 >        }
501 >        ii++;
502 >      }
503 >      
504 >      
505 >      if (foundNBI.size() > 0) {
506 >        // sort the foundNBI by the score:
507 >        std::sort(foundNBI.begin(), foundNBI.end());
508 >        
509 >        int bestScore = foundNBI[0].first;
510 >        std::vector<std::string> theKeys = foundNBI[0].second;
511 >        
512 >        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
513 >        return bestType;
514 >      } else {
515 >        //if no exact match found, try wild card match
516 >        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
517 >      }
518 >    }
519    }
520    
521    BondType* ForceField::getExactBondType(const std::string &at1,
# Line 390 | Line 547 | namespace oopse {
547      keys.push_back(at4);  
548      return torsionTypeCont_.find(keys);
549    }
550 <
550 >  
551 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
552 >                                                   const std::string &at2,
553 >                                                   const std::string &at3,
554 >                                                   const std::string &at4){
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_.find(keys);
561 >  }
562 >  
563    NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
564      std::vector<std::string> keys;
565      keys.push_back(at1);
566      keys.push_back(at2);    
567      return nonBondedInteractionTypeCont_.find(keys);
568    }
569 +  
570  
401
571    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
572      std::vector<std::string> keys;
573      keys.push_back(at);
574 +    atypeIdentToName[atomType->getIdent()] = at;
575      return atomTypeCont_.add(keys, atomType);
576    }
577  
578 +  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
579 +    std::vector<std::string> keys;
580 +    keys.push_back(at);
581 +    atypeIdentToName[atomType->getIdent()] = at;
582 +    return atomTypeCont_.replace(keys, atomType);
583 +  }
584 +
585    bool ForceField::addBondType(const std::string &at1, const std::string &at2,
586                                 BondType* bondType) {
587      std::vector<std::string> keys;
# Line 435 | Line 612 | namespace oopse {
612      return torsionTypeCont_.add(keys, torsionType);
613    }
614  
615 +  bool ForceField::addInversionType(const std::string &at1,
616 +                                    const std::string &at2,
617 +                                    const std::string &at3,
618 +                                    const std::string &at4,
619 +                                    InversionType* inversionType) {
620 +    std::vector<std::string> keys;
621 +    keys.push_back(at1);
622 +    keys.push_back(at2);    
623 +    keys.push_back(at3);    
624 +    keys.push_back(at4);    
625 +    return inversionTypeCont_.add(keys, inversionType);
626 +  }
627 +  
628    bool ForceField::addNonBondedInteractionType(const std::string &at1,
629                                                 const std::string &at2,
630                                                 NonBondedInteractionType* nbiType) {
# Line 463 | Line 653 | namespace oopse {
653          } else {
654            sprintf( painCave.errMsg,
655                     "Can not cast GenericData to LJParam\n");
656 <          painCave.severity = OOPSE_ERROR;
656 >          painCave.severity = OPENMD_ERROR;
657            painCave.isFatal = 1;
658            simError();          
659          }            
660        } else {
661          sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
662 <        painCave.severity = OOPSE_ERROR;
662 >        painCave.severity = OPENMD_ERROR;
663          painCave.isFatal = 1;
664          simError();          
665        }
# Line 499 | Line 689 | namespace oopse {
689                   "\tHave you tried setting the FORCE_PARAM_PATH environment "
690                   "variable?\n",
691                   forceFieldFilename.c_str() );
692 <        painCave.severity = OOPSE_ERROR;
692 >        painCave.severity = OPENMD_ERROR;
693          painCave.isFatal = 1;
694          simError();
695        }
# Line 507 | Line 697 | namespace oopse {
697      return ffStream;
698    }
699  
700 <  void ForceField::setFortranForceOptions(){
511 <    ForceOptions theseFortranOptions;
512 <    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
513 <    setfForceOptions(&theseFortranOptions);
514 <  }
515 < } //end namespace oopse
700 > } //end namespace OpenMD

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 1269 by gezelter, Tue Jul 1 13:28:23 2008 UTC vs.
branches/development/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines