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

Comparing:
branches/development/src/UseTheForce/ForceField.cpp (file contents), Revision 1473 by gezelter, Tue Jul 20 15:43:00 2010 UTC vs.
trunk/src/brains/ForceField.cpp (file contents), Revision 1880 by gezelter, Mon Jun 17 18:28:30 2013 UTC

# Line 35 | Line 35
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).                        
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"
53 >
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() {
83 >  ForceField::ForceField(std::string ffName) {
84  
85      char* tempPath;
86      tempPath = getenv("FORCE_PARAM_PATH");
# Line 67 | Line 91 | namespace OpenMD {
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 115 | Line 227 | namespace OpenMD {
227  
228        std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
229  
230 <      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
230 >      for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
231          jj = 0;
232 <        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
232 >        for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
233  
234            bondTypeScore = ii + jj;
235  
# Line 135 | Line 247 | namespace OpenMD {
247        }
248  
249  
250 <      if (foundBonds.size() > 0) {
250 >      if (!foundBonds.empty()) {
251          // sort the foundBonds by the score:
252          std::sort(foundBonds.begin(), foundBonds.end());
253      
142 –        int bestScore = foundBonds[0].first;
254          std::vector<std::string> theKeys = foundBonds[0].second;
255          
256          BondType* bestType = bondTypeCont_.find(theKeys);
# Line 197 | Line 308 | namespace OpenMD {
308  
309        std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
310  
311 <      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
311 >      for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
312          ii = 0;
313 <        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
313 >        for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
314            kk = 0;
315 <          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
315 >          for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
316            
317              IKscore = ii + kk;
318  
# Line 221 | Line 332 | namespace OpenMD {
332          jj++;
333        }
334        
335 <      if (foundBends.size() > 0) {
335 >      if (!foundBends.empty()) {
336          std::sort(foundBends.begin(), foundBends.end());
226 –        int jscore = foundBends[0].first;
227 –        int ikscore = foundBends[0].second;
337          std::vector<std::string> theKeys = foundBends[0].third;      
338          
339          BendType* bestType = bendTypeCont_.find(theKeys);  
# Line 293 | Line 402 | namespace OpenMD {
402  
403        std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
404  
405 <      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
405 >      for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
406          kk = 0;
407 <        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
407 >        for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
408            ii = 0;      
409 <          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
409 >          for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
410              ll = 0;
411 <            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
411 >            for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
412            
413                ILscore = ii + ll;
414                JKscore = jj + kk;
# Line 323 | Line 432 | namespace OpenMD {
432          jj++;
433        }
434        
435 <      if (foundTorsions.size() > 0) {
435 >      if (!foundTorsions.empty()) {
436          std::sort(foundTorsions.begin(), foundTorsions.end());
328 –        int jkscore = foundTorsions[0].first;
329 –        int ilscore = foundTorsions[0].second;
437          std::vector<std::string> theKeys = foundTorsions[0].third;
438          
439          TorsionType* bestType = torsionTypeCont_.find(theKeys);
# Line 394 | Line 501 | namespace OpenMD {
501        
502        std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
503        
504 <      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
504 >      for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
505          kk = 0;
506 <        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
506 >        for (k = at3Chain.begin(); k != at3Chain.end(); ++k) {
507            ii = 0;      
508 <          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
508 >          for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
509              ll = 0;
510 <            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
510 >            for (l = at4Chain.begin(); l != at4Chain.end(); ++l) {
511                
512                Iscore = ii;
513                JKLscore = jj + kk + ll;
# Line 424 | Line 531 | namespace OpenMD {
531          jj++;
532        }
533          
534 <      if (foundInversions.size() > 0) {
534 >      if (!foundInversions.empty()) {
535          std::sort(foundInversions.begin(), foundInversions.end());
429 –        int iscore = foundInversions[0].first;
430 –        int jklscore = foundInversions[0].second;
536          std::vector<std::string> theKeys = foundInversions[0].third;
537          
538          InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
# Line 440 | Line 545 | namespace OpenMD {
545    }
546    
547    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
548 +    
549      std::vector<std::string> keys;
550      keys.push_back(at1);
551      keys.push_back(at2);    
# Line 449 | Line 555 | namespace OpenMD {
555      if (nbiType) {
556        return nbiType;
557      } else {
558 <      //if no exact match found, try wild card match
559 <      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
560 <    }    
558 >      AtomType* atype1;
559 >      AtomType* atype2;
560 >      std::vector<std::string> at1key;
561 >      at1key.push_back(at1);
562 >      atype1 = atomTypeCont_.find(at1key);
563 >      
564 >      std::vector<std::string> at2key;
565 >      at2key.push_back(at2);
566 >      atype2 = atomTypeCont_.find(at2key);
567 >      
568 >      // query atom types for their chains of responsibility
569 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
570 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
571 >      
572 >      std::vector<AtomType*>::iterator i;
573 >      std::vector<AtomType*>::iterator j;
574 >      
575 >      int ii = 0;
576 >      int jj = 0;
577 >      int nbiTypeScore;
578 >      
579 >      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
580 >      
581 >      for (i = at1Chain.begin(); i != at1Chain.end(); ++i) {
582 >        jj = 0;
583 >        for (j = at2Chain.begin(); j != at2Chain.end(); ++j) {
584 >          
585 >          nbiTypeScore = ii + jj;
586 >          
587 >          std::vector<std::string> myKeys;
588 >          myKeys.push_back((*i)->getName());
589 >          myKeys.push_back((*j)->getName());
590 >          
591 >          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
592 >          if (nbiType) {
593 >            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
594 >          }
595 >          jj++;
596 >        }
597 >        ii++;
598 >      }
599 >      
600 >      
601 >      if (!foundNBI.empty()) {
602 >        // sort the foundNBI by the score:
603 >        std::sort(foundNBI.begin(), foundNBI.end());      
604 >        std::vector<std::string> theKeys = foundNBI[0].second;
605 >        
606 >        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
607 >        return bestType;
608 >      } else {
609 >        //if no exact match found, try wild card match
610 >        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
611 >      }
612 >    }
613    }
614    
615    BondType* ForceField::getExactBondType(const std::string &at1,
# Line 507 | Line 665 | namespace OpenMD {
665    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
666      std::vector<std::string> keys;
667      keys.push_back(at);
668 +    atypeIdentToName[atomType->getIdent()] = at;
669      return atomTypeCont_.add(keys, atomType);
670    }
671  
672    bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
673      std::vector<std::string> keys;
674      keys.push_back(at);
675 +    atypeIdentToName[atomType->getIdent()] = at;
676      return atomTypeCont_.replace(keys, atomType);
677    }
678  
# Line 569 | Line 729 | namespace OpenMD {
729    }
730    
731    RealType ForceField::getRcutFromAtomType(AtomType* at) {
732 <    /**@todo */
573 <    GenericData* data;
574 <    RealType rcut = 0.0;
732 >    RealType rcut(0.0);
733      
734 <    if (at->isLennardJones()) {
735 <      data = at->getPropertyByName("LennardJones");
736 <      if (data != NULL) {
579 <        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
580 <        
581 <        if (ljData != NULL) {
582 <          LJParam ljParam = ljData->getData();
583 <          
584 <          //by default use 2.5*sigma as cutoff radius
585 <          rcut = 2.5 * ljParam.sigma;
586 <          
587 <        } else {
588 <          sprintf( painCave.errMsg,
589 <                   "Can not cast GenericData to LJParam\n");
590 <          painCave.severity = OPENMD_ERROR;
591 <          painCave.isFatal = 1;
592 <          simError();          
593 <        }            
594 <      } else {
595 <        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
596 <        painCave.severity = OPENMD_ERROR;
597 <        painCave.isFatal = 1;
598 <        simError();          
599 <      }
734 >    LennardJonesAdapter lja = LennardJonesAdapter(at);
735 >    if (lja.isLennardJones()) {
736 >      rcut = 2.5 * lja.getSigma();
737      }
738 +    EAMAdapter ea = EAMAdapter(at);
739 +    if (ea.isEAM()) {
740 +      rcut = max(rcut, ea.getRcut());
741 +    }
742 +    SuttonChenAdapter sca = SuttonChenAdapter(at);
743 +    if (sca.isSuttonChen()) {
744 +      rcut = max(rcut, 2.0 * sca.getAlpha());
745 +    }
746 +    GayBerneAdapter gba = GayBerneAdapter(at);
747 +    if (gba.isGayBerne()) {
748 +      rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
749 +    }
750 +    StickyAdapter sa = StickyAdapter(at);
751 +    if (sa.isSticky()) {
752 +      rcut = max(rcut, max(sa.getRu(), sa.getRup()));
753 +    }
754 +
755      return rcut;    
756    }
757    
# Line 631 | Line 785 | namespace OpenMD {
785      return ffStream;
786    }
787  
634 –  void ForceField::setFortranForceOptions(){
635 –    ForceOptions theseFortranOptions;
636 –    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
637 –    setfForceOptions(&theseFortranOptions);
638 –  }
788   } //end namespace OpenMD

Diff Legend

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