ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/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/brains/ForceField.cpp (file contents), Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 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, 234107 (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
46 * @time 22:51am
47   * @version 1.0
48   */
49    
50   #include <algorithm>
51 < #include "UseTheForce/ForceField.hpp"
51 > #include "brains/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 {
53  
54 <  ForceField::ForceField() {
54 > #include "io/OptionSectionParser.hpp"
55 > #include "io/BaseAtomTypesSectionParser.hpp"
56 > #include "io/DirectionalAtomTypesSectionParser.hpp"
57 > #include "io/AtomTypesSectionParser.hpp"
58 > #include "io/BendTypesSectionParser.hpp"
59 > #include "io/BondTypesSectionParser.hpp"
60 > #include "io/ChargeAtomTypesSectionParser.hpp"
61 > #include "io/EAMAtomTypesSectionParser.hpp"
62 > #include "io/FluctuatingChargeAtomTypesSectionParser.hpp"
63 > #include "io/GayBerneAtomTypesSectionParser.hpp"
64 > #include "io/InversionTypesSectionParser.hpp"
65 > #include "io/LennardJonesAtomTypesSectionParser.hpp"
66 > #include "io/MultipoleAtomTypesSectionParser.hpp"
67 > #include "io/NonBondedInteractionsSectionParser.hpp"
68 > #include "io/PolarizableAtomTypesSectionParser.hpp"
69 > #include "io/SCAtomTypesSectionParser.hpp"
70 > #include "io/ShapeAtomTypesSectionParser.hpp"
71 > #include "io/StickyAtomTypesSectionParser.hpp"
72 > #include "io/StickyPowerAtomTypesSectionParser.hpp"
73 > #include "io/TorsionTypesSectionParser.hpp"
74 >
75 > #include "types/LennardJonesAdapter.hpp"
76 > #include "types/EAMAdapter.hpp"
77 > #include "types/SuttonChenAdapter.hpp"
78 > #include "types/GayBerneAdapter.hpp"
79 > #include "types/StickyAdapter.hpp"
80 >
81 > namespace OpenMD {
82 >
83 >  ForceField::ForceField(std::string ffName) {
84 >
85      char* tempPath;
86      tempPath = getenv("FORCE_PARAM_PATH");
87 <
87 >    
88      if (tempPath == NULL) {
89        //convert a macro from compiler to a string in c++
90        STR_DEFINE(ffPath_, FRC_PATH );
91      } else {
92        ffPath_ = tempPath;
93      }
94 +
95 +    setForceFieldFileName(ffName + ".frc");
96 +
97 +    /**
98 +     * The order of adding section parsers is important.
99 +     *
100 +     * OptionSectionParser must come first to set options for other
101 +     * parsers
102 +     *
103 +     * DirectionalAtomTypesSectionParser should be added before
104 +     * AtomTypesSectionParser, and these two section parsers will
105 +     * actually create "real" AtomTypes (AtomTypesSectionParser will
106 +     * create AtomType and DirectionalAtomTypesSectionParser will
107 +     * create DirectionalAtomType, which is a subclass of AtomType and
108 +     * should come first).
109 +     *
110 +     * Other AtomTypes Section Parsers will not create the "real"
111 +     * AtomType, they only add and set some attributes of the AtomType
112 +     * (via the Adapters). Thus ordering of these is not important.
113 +     * AtomTypesSectionParser should be added before other atom type
114 +     *
115 +     * The order of BondTypesSectionParser, BendTypesSectionParser and
116 +     * TorsionTypesSectionParser, etc. are not important.
117 +     */
118 +
119 +    spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
120 +    spMan_.push_back(new BaseAtomTypesSectionParser());
121 +    spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
122 +    spMan_.push_back(new AtomTypesSectionParser());
123 +
124 +    spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
125 +    spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
126 +    spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
127 +    spMan_.push_back(new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_));
128 +    spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_));
129 +    spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
130 +    spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
131 +    spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
132 +    spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_));
133 +    spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
134 +    spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
135 +
136 +    spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
137 +    spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
138 +    spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
139 +    spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_));
140 +
141 +    spMan_.push_back(new NonBondedInteractionsSectionParser(forceFieldOptions_));    
142    }
143  
144 +  void ForceField::parse(const std::string& filename) {
145 +    ifstrstream* ffStream;
146  
147 <  ForceField::~ForceField() {
148 <    deleteAtypes();
149 <    deleteSwitch();
147 >    ffStream = openForceFieldFile(filename);
148 >
149 >    spMan_.parse(*ffStream, *this);
150 >
151 >    ForceField::AtomTypeContainer::MapTypeIterator i;
152 >    AtomType* at;
153 >
154 >    for (at = atomTypeCont_.beginType(i); at != NULL;
155 >         at = atomTypeCont_.nextType(i)) {
156 >
157 >      // useBase sets the responsibilities, and these have to be done
158 >      // after the atomTypes and Base types have all been scanned:
159 >
160 >      std::vector<AtomType*> ayb = at->allYourBase();      
161 >      if (ayb.size() > 1) {
162 >        for (int j = ayb.size()-1; j > 0; j--) {
163 >          
164 >          ayb[j-1]->useBase(ayb[j]);
165 >
166 >        }
167 >      }
168 >    }
169 >
170 >    delete ffStream;
171    }
172  
173 +  /**
174 +   * getAtomType by string
175 +   *
176 +   * finds the requested atom type in this force field using the string
177 +   * name of the atom type.
178 +   */
179    AtomType* ForceField::getAtomType(const std::string &at) {
180      std::vector<std::string> keys;
181      keys.push_back(at);
182      return atomTypeCont_.find(keys);
183    }
184  
185 +  /**
186 +   * getAtomType by ident
187 +   *
188 +   * finds the requested atom type in this force field using the
189 +   * integer ident instead of the string name of the atom type.
190 +   */
191 +  AtomType* ForceField::getAtomType(int ident) {  
192 +    std::string at = atypeIdentToName.find(ident)->second;
193 +    return getAtomType(at);
194 +  }
195 +
196    BondType* ForceField::getBondType(const std::string &at1,
197                                      const std::string &at2) {
198      std::vector<std::string> keys;
# Line 133 | Line 246 | namespace oopse {
246          ii++;
247        }
248  
136      // sort the foundBonds by the score:
249  
250 <      std::sort(foundBonds.begin(), foundBonds.end());
251 <      
252 <      int bestScore = foundBonds[0].first;
253 <      std::vector<std::string> theKeys = foundBonds[0].second;
254 <
255 <      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
256 <      BondType* bestType = bondTypeCont_.find(theKeys);
257 <      if (bestType)
258 <        return bestType;
259 <      else {
260 <        //if no exact match found, try wild card match
261 <        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
250 >      if (foundBonds.size() > 0) {
251 >        // sort the foundBonds by the score:
252 >        std::sort(foundBonds.begin(), foundBonds.end());
253 >    
254 >        int bestScore = foundBonds[0].first;
255 >        std::vector<std::string> theKeys = foundBonds[0].second;
256 >        
257 >        BondType* bestType = bondTypeCont_.find(theKeys);
258 >        
259 >        return bestType;
260 >      } else {
261 >        //if no exact match found, try wild card match
262 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
263        }
264      }
265    }
# Line 220 | Line 333 | namespace oopse {
333          jj++;
334        }
335        
336 <      std::sort(foundBends.begin(), foundBends.end());
337 <
338 <      int jscore = foundBends[0].first;
339 <      int ikscore = foundBends[0].second;
340 <      std::vector<std::string> theKeys = foundBends[0].third;
341 <
342 <      std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
343 <
344 <      BendType* bestType = bendTypeCont_.find(theKeys);  
232 <      if (bestType)
233 <        return bestType;
234 <      else {
235 <      
336 >      if (foundBends.size() > 0) {
337 >        std::sort(foundBends.begin(), foundBends.end());
338 >        int jscore = foundBends[0].first;
339 >        int ikscore = foundBends[0].second;
340 >        std::vector<std::string> theKeys = foundBends[0].third;      
341 >        
342 >        BendType* bestType = bendTypeCont_.find(theKeys);  
343 >        return bestType;
344 >      } else {        
345          //if no exact match found, try wild card match
346          return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
347        }
348      }
349    }
350  
242
351    TorsionType* ForceField::getTorsionType(const std::string &at1,
352                                            const std::string &at2,
353                                            const std::string &at3,
# Line 327 | Line 435 | namespace oopse {
435          jj++;
436        }
437        
438 <      std::sort(foundTorsions.begin(), foundTorsions.end());
438 >      if (foundTorsions.size() > 0) {
439 >        std::sort(foundTorsions.begin(), foundTorsions.end());
440 >        int jkscore = foundTorsions[0].first;
441 >        int ilscore = foundTorsions[0].second;
442 >        std::vector<std::string> theKeys = foundTorsions[0].third;
443 >        
444 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
445 >        return bestType;
446 >      } else {
447 >        //if no exact match found, try wild card match
448 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
449 >      }
450 >    }
451 >  }
452  
453 <      int jkscore = foundTorsions[0].first;
454 <      int ilscore = foundTorsions[0].second;
455 <      std::vector<std::string> theKeys = foundTorsions[0].third;
453 >  InversionType* ForceField::getInversionType(const std::string &at1,
454 >                                              const std::string &at2,
455 >                                              const std::string &at3,
456 >                                              const std::string &at4) {
457 >    std::vector<std::string> keys;
458 >    keys.push_back(at1);
459 >    keys.push_back(at2);    
460 >    keys.push_back(at3);    
461 >    keys.push_back(at4);    
462  
463 <      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
463 >    //try exact match first
464 >    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
465 >    if (inversionType) {
466 >      return inversionType;
467 >    } else {
468 >      
469 >      AtomType* atype1;
470 >      AtomType* atype2;
471 >      AtomType* atype3;
472 >      AtomType* atype4;
473 >      std::vector<std::string> at1key;
474 >      at1key.push_back(at1);
475 >      atype1 = atomTypeCont_.find(at1key);
476 >      
477 >      std::vector<std::string> at2key;
478 >      at2key.push_back(at2);
479 >      atype2 = atomTypeCont_.find(at2key);
480 >      
481 >      std::vector<std::string> at3key;
482 >      at3key.push_back(at3);
483 >      atype3 = atomTypeCont_.find(at3key);
484 >      
485 >      std::vector<std::string> at4key;
486 >      at4key.push_back(at4);
487 >      atype4 = atomTypeCont_.find(at4key);
488  
489 +      // query atom types for their chains of responsibility
490 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
491 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
492 +      std::vector<AtomType*> at3Chain = atype3->allYourBase();
493 +      std::vector<AtomType*> at4Chain = atype4->allYourBase();
494 +
495 +      std::vector<AtomType*>::iterator i;
496 +      std::vector<AtomType*>::iterator j;
497 +      std::vector<AtomType*>::iterator k;
498 +      std::vector<AtomType*>::iterator l;
499 +
500 +      int ii = 0;
501 +      int jj = 0;
502 +      int kk = 0;
503 +      int ll = 0;
504 +      int Iscore;
505 +      int JKLscore;
506        
507 <      TorsionType* bestType = torsionTypeCont_.find(theKeys);
508 <      if (bestType) {
509 <        return bestType;
507 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
508 >      
509 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
510 >        kk = 0;
511 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
512 >          ii = 0;      
513 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
514 >            ll = 0;
515 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
516 >              
517 >              Iscore = ii;
518 >              JKLscore = jj + kk + ll;
519 >              
520 >              std::vector<std::string> myKeys;
521 >              myKeys.push_back((*i)->getName());
522 >              myKeys.push_back((*j)->getName());
523 >              myKeys.push_back((*k)->getName());
524 >              myKeys.push_back((*l)->getName());
525 >              
526 >              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
527 >              if (inversionType) {
528 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
529 >              }
530 >              ll++;
531 >            }
532 >            ii++;
533 >          }
534 >          kk++;
535 >        }
536 >        jj++;
537 >      }
538 >        
539 >      if (foundInversions.size() > 0) {
540 >        std::sort(foundInversions.begin(), foundInversions.end());
541 >        int iscore = foundInversions[0].first;
542 >        int jklscore = foundInversions[0].second;
543 >        std::vector<std::string> theKeys = foundInversions[0].third;
544 >        
545 >        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
546 >        return bestType;
547        } else {
548          //if no exact match found, try wild card match
549 <        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
549 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
550        }
551      }
552    }
553 <
553 >  
554    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
555 +    
556      std::vector<std::string> keys;
557      keys.push_back(at1);
558      keys.push_back(at2);    
# Line 356 | Line 562 | namespace oopse {
562      if (nbiType) {
563        return nbiType;
564      } else {
565 <      //if no exact match found, try wild card match
566 <      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
567 <    }    
565 >      AtomType* atype1;
566 >      AtomType* atype2;
567 >      std::vector<std::string> at1key;
568 >      at1key.push_back(at1);
569 >      atype1 = atomTypeCont_.find(at1key);
570 >      
571 >      std::vector<std::string> at2key;
572 >      at2key.push_back(at2);
573 >      atype2 = atomTypeCont_.find(at2key);
574 >      
575 >      // query atom types for their chains of responsibility
576 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
577 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
578 >      
579 >      std::vector<AtomType*>::iterator i;
580 >      std::vector<AtomType*>::iterator j;
581 >      
582 >      int ii = 0;
583 >      int jj = 0;
584 >      int nbiTypeScore;
585 >      
586 >      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
587 >      
588 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
589 >        jj = 0;
590 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
591 >          
592 >          nbiTypeScore = ii + jj;
593 >          
594 >          std::vector<std::string> myKeys;
595 >          myKeys.push_back((*i)->getName());
596 >          myKeys.push_back((*j)->getName());
597 >          
598 >          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
599 >          if (nbiType) {
600 >            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
601 >          }
602 >          jj++;
603 >        }
604 >        ii++;
605 >      }
606 >      
607 >      
608 >      if (foundNBI.size() > 0) {
609 >        // sort the foundNBI by the score:
610 >        std::sort(foundNBI.begin(), foundNBI.end());
611 >        
612 >        int bestScore = foundNBI[0].first;
613 >        std::vector<std::string> theKeys = foundNBI[0].second;
614 >        
615 >        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
616 >        return bestType;
617 >      } else {
618 >        //if no exact match found, try wild card match
619 >        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
620 >      }
621 >    }
622    }
623    
624    BondType* ForceField::getExactBondType(const std::string &at1,
# Line 390 | Line 650 | namespace oopse {
650      keys.push_back(at4);  
651      return torsionTypeCont_.find(keys);
652    }
653 <
653 >  
654 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
655 >                                                   const std::string &at2,
656 >                                                   const std::string &at3,
657 >                                                   const std::string &at4){
658 >    std::vector<std::string> keys;
659 >    keys.push_back(at1);
660 >    keys.push_back(at2);    
661 >    keys.push_back(at3);    
662 >    keys.push_back(at4);  
663 >    return inversionTypeCont_.find(keys);
664 >  }
665 >  
666    NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
667      std::vector<std::string> keys;
668      keys.push_back(at1);
669      keys.push_back(at2);    
670      return nonBondedInteractionTypeCont_.find(keys);
671    }
672 +  
673  
401
674    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
675      std::vector<std::string> keys;
676      keys.push_back(at);
677 +    atypeIdentToName[atomType->getIdent()] = at;
678      return atomTypeCont_.add(keys, atomType);
679    }
680  
681 +  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
682 +    std::vector<std::string> keys;
683 +    keys.push_back(at);
684 +    atypeIdentToName[atomType->getIdent()] = at;
685 +    return atomTypeCont_.replace(keys, atomType);
686 +  }
687 +
688    bool ForceField::addBondType(const std::string &at1, const std::string &at2,
689                                 BondType* bondType) {
690      std::vector<std::string> keys;
# Line 435 | Line 715 | namespace oopse {
715      return torsionTypeCont_.add(keys, torsionType);
716    }
717  
718 +  bool ForceField::addInversionType(const std::string &at1,
719 +                                    const std::string &at2,
720 +                                    const std::string &at3,
721 +                                    const std::string &at4,
722 +                                    InversionType* inversionType) {
723 +    std::vector<std::string> keys;
724 +    keys.push_back(at1);
725 +    keys.push_back(at2);    
726 +    keys.push_back(at3);    
727 +    keys.push_back(at4);    
728 +    return inversionTypeCont_.add(keys, inversionType);
729 +  }
730 +  
731    bool ForceField::addNonBondedInteractionType(const std::string &at1,
732                                                 const std::string &at2,
733                                                 NonBondedInteractionType* nbiType) {
# Line 445 | Line 738 | namespace oopse {
738    }
739    
740    RealType ForceField::getRcutFromAtomType(AtomType* at) {
741 <    /**@todo */
449 <    GenericData* data;
450 <    RealType rcut = 0.0;
741 >    RealType rcut(0.0);
742      
743 <    if (at->isLennardJones()) {
744 <      data = at->getPropertyByName("LennardJones");
745 <      if (data != NULL) {
455 <        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
456 <        
457 <        if (ljData != NULL) {
458 <          LJParam ljParam = ljData->getData();
459 <          
460 <          //by default use 2.5*sigma as cutoff radius
461 <          rcut = 2.5 * ljParam.sigma;
462 <          
463 <        } else {
464 <          sprintf( painCave.errMsg,
465 <                   "Can not cast GenericData to LJParam\n");
466 <          painCave.severity = OOPSE_ERROR;
467 <          painCave.isFatal = 1;
468 <          simError();          
469 <        }            
470 <      } else {
471 <        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
472 <        painCave.severity = OOPSE_ERROR;
473 <        painCave.isFatal = 1;
474 <        simError();          
475 <      }
743 >    LennardJonesAdapter lja = LennardJonesAdapter(at);
744 >    if (lja.isLennardJones()) {
745 >      rcut = 2.5 * lja.getSigma();
746      }
747 +    EAMAdapter ea = EAMAdapter(at);
748 +    if (ea.isEAM()) {
749 +      rcut = max(rcut, ea.getRcut());
750 +    }
751 +    SuttonChenAdapter sca = SuttonChenAdapter(at);
752 +    if (sca.isSuttonChen()) {
753 +      rcut = max(rcut, 2.0 * sca.getAlpha());
754 +    }
755 +    GayBerneAdapter gba = GayBerneAdapter(at);
756 +    if (gba.isGayBerne()) {
757 +      rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
758 +    }
759 +    StickyAdapter sa = StickyAdapter(at);
760 +    if (sa.isSticky()) {
761 +      rcut = max(rcut, max(sa.getRu(), sa.getRup()));
762 +    }
763 +
764      return rcut;    
765    }
766    
# Line 499 | Line 786 | namespace oopse {
786                   "\tHave you tried setting the FORCE_PARAM_PATH environment "
787                   "variable?\n",
788                   forceFieldFilename.c_str() );
789 <        painCave.severity = OOPSE_ERROR;
789 >        painCave.severity = OPENMD_ERROR;
790          painCave.isFatal = 1;
791          simError();
792        }
# Line 507 | Line 794 | namespace oopse {
794      return ffStream;
795    }
796  
797 <  void ForceField::setFortranForceOptions(){
511 <    ForceOptions theseFortranOptions;
512 <    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
513 <    setfForceOptions(&theseFortranOptions);
514 <  }
515 < } //end namespace oopse
797 > } //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/brains/ForceField.cpp (property svn:keywords), Revision 1850 by gezelter, Wed Feb 20 15:39:39 2013 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines