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 1530 by gezelter, Tue Dec 28 21:47:55 2010 UTC vs.
trunk/src/brains/ForceField.cpp (file contents), Revision 1790 by gezelter, Thu Aug 30 17:18:22 2012 UTC

# Line 36 | Line 36
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).                        
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 48 | Line 49
49   */
50    
51   #include <algorithm>
52 < #include "UseTheForce/ForceField.hpp"
52 > #include "brains/ForceField.hpp"
53   #include "utils/simError.h"
54 < #include "utils/Tuple.hpp"
55 < #include "UseTheForce/DarkSide/atype_interface.h"
54 >
55 > #include "io/OptionSectionParser.hpp"
56 > #include "io/BaseAtomTypesSectionParser.hpp"
57 > #include "io/DirectionalAtomTypesSectionParser.hpp"
58 > #include "io/AtomTypesSectionParser.hpp"
59 > #include "io/BendTypesSectionParser.hpp"
60 > #include "io/BondTypesSectionParser.hpp"
61 > #include "io/ChargeAtomTypesSectionParser.hpp"
62 > #include "io/EAMAtomTypesSectionParser.hpp"
63 > #include "io/FluctuatingChargeAtomTypesSectionParser.hpp"
64 > #include "io/GayBerneAtomTypesSectionParser.hpp"
65 > #include "io/InversionTypesSectionParser.hpp"
66 > #include "io/LennardJonesAtomTypesSectionParser.hpp"
67 > #include "io/MultipoleAtomTypesSectionParser.hpp"
68 > #include "io/NonBondedInteractionsSectionParser.hpp"
69 > #include "io/PolarizableAtomTypesSectionParser.hpp"
70 > #include "io/SCAtomTypesSectionParser.hpp"
71 > #include "io/ShapeAtomTypesSectionParser.hpp"
72 > #include "io/StickyAtomTypesSectionParser.hpp"
73 > #include "io/StickyPowerAtomTypesSectionParser.hpp"
74 > #include "io/TorsionTypesSectionParser.hpp"
75 >
76 > #include "types/LennardJonesAdapter.hpp"
77 > #include "types/EAMAdapter.hpp"
78 > #include "types/SuttonChenAdapter.hpp"
79 > #include "types/GayBerneAdapter.hpp"
80 > #include "types/StickyAdapter.hpp"
81 >
82   namespace OpenMD {
83  
84 <  ForceField::ForceField() {
84 >  ForceField::ForceField(std::string ffName) {
85  
86      char* tempPath;
87      tempPath = getenv("FORCE_PARAM_PATH");
# Line 65 | Line 92 | namespace OpenMD {
92      } else {
93        ffPath_ = tempPath;
94      }
95 +
96 +    setForceFieldFileName(ffName + ".frc");
97 +
98 +    /**
99 +     * The order of adding section parsers is important.
100 +     *
101 +     * OptionSectionParser must come first to set options for other
102 +     * parsers
103 +     *
104 +     * DirectionalAtomTypesSectionParser should be added before
105 +     * AtomTypesSectionParser, and these two section parsers will
106 +     * actually create "real" AtomTypes (AtomTypesSectionParser will
107 +     * create AtomType and DirectionalAtomTypesSectionParser will
108 +     * create DirectionalAtomType, which is a subclass of AtomType and
109 +     * should come first).
110 +     *
111 +     * Other AtomTypes Section Parsers will not create the "real"
112 +     * AtomType, they only add and set some attributes of the AtomType
113 +     * (via the Adapters). Thus ordering of these is not important.
114 +     * AtomTypesSectionParser should be added before other atom type
115 +     *
116 +     * The order of BondTypesSectionParser, BendTypesSectionParser and
117 +     * TorsionTypesSectionParser, etc. are not important.
118 +     */
119 +
120 +    spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
121 +    spMan_.push_back(new BaseAtomTypesSectionParser());
122 +    spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
123 +    spMan_.push_back(new AtomTypesSectionParser());
124 +
125 +    spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
126 +    spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
127 +    spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
128 +    spMan_.push_back(new FluctuatingChargeAtomTypesSectionParser(forceFieldOptions_));
129 +    spMan_.push_back(new PolarizableAtomTypesSectionParser(forceFieldOptions_));
130 +    spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
131 +    spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
132 +    spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
133 +    spMan_.push_back(new ShapeAtomTypesSectionParser(forceFieldOptions_));
134 +    spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
135 +    spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
136 +
137 +    spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
138 +    spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
139 +    spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
140 +    spMan_.push_back(new InversionTypesSectionParser(forceFieldOptions_));
141 +
142 +    spMan_.push_back(new NonBondedInteractionsSectionParser(forceFieldOptions_));    
143    }
144  
145 +  void ForceField::parse(const std::string& filename) {
146 +    ifstrstream* ffStream;
147  
148 <  ForceField::~ForceField() {
149 <    deleteAtypes();
148 >    ffStream = openForceFieldFile(filename);
149 >
150 >    spMan_.parse(*ffStream, *this);
151 >
152 >    ForceField::AtomTypeContainer::MapTypeIterator i;
153 >    AtomType* at;
154 >
155 >    for (at = atomTypeCont_.beginType(i); at != NULL;
156 >         at = atomTypeCont_.nextType(i)) {
157 >
158 >      // useBase sets the responsibilities, and these have to be done
159 >      // after the atomTypes and Base types have all been scanned:
160 >
161 >      std::vector<AtomType*> ayb = at->allYourBase();      
162 >      if (ayb.size() > 1) {
163 >        for (int j = ayb.size()-1; j > 0; j--) {
164 >          
165 >          ayb[j-1]->useBase(ayb[j]);
166 >
167 >        }
168 >      }
169 >    }
170 >
171 >    delete ffStream;
172    }
173  
174 +  /**
175 +   * getAtomType by string
176 +   *
177 +   * finds the requested atom type in this force field using the string
178 +   * name of the atom type.
179 +   */
180    AtomType* ForceField::getAtomType(const std::string &at) {
181      std::vector<std::string> keys;
182      keys.push_back(at);
183      return atomTypeCont_.find(keys);
184    }
185  
186 +  /**
187 +   * getAtomType by ident
188 +   *
189 +   * finds the requested atom type in this force field using the
190 +   * integer ident instead of the string name of the atom type.
191 +   */
192 +  AtomType* ForceField::getAtomType(int ident) {  
193 +    std::string at = atypeIdentToName.find(ident)->second;
194 +    return getAtomType(at);
195 +  }
196 +
197    BondType* ForceField::getBondType(const std::string &at1,
198                                      const std::string &at2) {
199      std::vector<std::string> keys;
# Line 437 | Line 553 | namespace OpenMD {
553    }
554    
555    NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
556 +    
557      std::vector<std::string> keys;
558      keys.push_back(at1);
559      keys.push_back(at2);    
# Line 446 | Line 563 | namespace OpenMD {
563      if (nbiType) {
564        return nbiType;
565      } else {
566 <      //if no exact match found, try wild card match
567 <      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
568 <    }    
566 >      AtomType* atype1;
567 >      AtomType* atype2;
568 >      std::vector<std::string> at1key;
569 >      at1key.push_back(at1);
570 >      atype1 = atomTypeCont_.find(at1key);
571 >      
572 >      std::vector<std::string> at2key;
573 >      at2key.push_back(at2);
574 >      atype2 = atomTypeCont_.find(at2key);
575 >      
576 >      // query atom types for their chains of responsibility
577 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
578 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
579 >      
580 >      std::vector<AtomType*>::iterator i;
581 >      std::vector<AtomType*>::iterator j;
582 >      
583 >      int ii = 0;
584 >      int jj = 0;
585 >      int nbiTypeScore;
586 >      
587 >      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
588 >      
589 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
590 >        jj = 0;
591 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
592 >          
593 >          nbiTypeScore = ii + jj;
594 >          
595 >          std::vector<std::string> myKeys;
596 >          myKeys.push_back((*i)->getName());
597 >          myKeys.push_back((*j)->getName());
598 >          
599 >          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
600 >          if (nbiType) {
601 >            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
602 >          }
603 >          jj++;
604 >        }
605 >        ii++;
606 >      }
607 >      
608 >      
609 >      if (foundNBI.size() > 0) {
610 >        // sort the foundNBI by the score:
611 >        std::sort(foundNBI.begin(), foundNBI.end());
612 >        
613 >        int bestScore = foundNBI[0].first;
614 >        std::vector<std::string> theKeys = foundNBI[0].second;
615 >        
616 >        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
617 >        return bestType;
618 >      } else {
619 >        //if no exact match found, try wild card match
620 >        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
621 >      }
622 >    }
623    }
624    
625    BondType* ForceField::getExactBondType(const std::string &at1,
# Line 504 | Line 675 | namespace OpenMD {
675    bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
676      std::vector<std::string> keys;
677      keys.push_back(at);
678 +    atypeIdentToName[atomType->getIdent()] = at;
679      return atomTypeCont_.add(keys, atomType);
680    }
681  
682    bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
683      std::vector<std::string> keys;
684      keys.push_back(at);
685 +    atypeIdentToName[atomType->getIdent()] = at;
686      return atomTypeCont_.replace(keys, atomType);
687    }
688  
# Line 566 | Line 739 | namespace OpenMD {
739    }
740    
741    RealType ForceField::getRcutFromAtomType(AtomType* at) {
742 <    /**@todo */
570 <    GenericData* data;
571 <    RealType rcut = 0.0;
742 >    RealType rcut(0.0);
743      
744 <    if (at->isLennardJones()) {
745 <      data = at->getPropertyByName("LennardJones");
746 <      if (data != NULL) {
576 <        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
577 <        
578 <        if (ljData != NULL) {
579 <          LJParam ljParam = ljData->getData();
580 <          
581 <          //by default use 2.5*sigma as cutoff radius
582 <          rcut = 2.5 * ljParam.sigma;
583 <          
584 <        } else {
585 <          sprintf( painCave.errMsg,
586 <                   "Can not cast GenericData to LJParam\n");
587 <          painCave.severity = OPENMD_ERROR;
588 <          painCave.isFatal = 1;
589 <          simError();          
590 <        }            
591 <      } else {
592 <        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
593 <        painCave.severity = OPENMD_ERROR;
594 <        painCave.isFatal = 1;
595 <        simError();          
596 <      }
744 >    LennardJonesAdapter lja = LennardJonesAdapter(at);
745 >    if (lja.isLennardJones()) {
746 >      rcut = 2.5 * lja.getSigma();
747      }
748 +    EAMAdapter ea = EAMAdapter(at);
749 +    if (ea.isEAM()) {
750 +      rcut = max(rcut, ea.getRcut());
751 +    }
752 +    SuttonChenAdapter sca = SuttonChenAdapter(at);
753 +    if (sca.isSuttonChen()) {
754 +      rcut = max(rcut, 2.0 * sca.getAlpha());
755 +    }
756 +    GayBerneAdapter gba = GayBerneAdapter(at);
757 +    if (gba.isGayBerne()) {
758 +      rcut = max(rcut, 2.5 * sqrt(2.0) * max(gba.getD(), gba.getL()));
759 +    }
760 +    StickyAdapter sa = StickyAdapter(at);
761 +    if (sa.isSticky()) {
762 +      rcut = max(rcut, max(sa.getRu(), sa.getRup()));
763 +    }
764 +
765      return rcut;    
766    }
767    
# Line 604 | Line 771 | namespace OpenMD {
771      ifstrstream* ffStream = new ifstrstream();
772      
773      //try to open the force filed file in current directory first    
774 <    ffStream->open(forceFieldFilename.c_str());
774 >    ffStream->open(forceFieldFilename.c_str(), ifstream::in | ifstream::binary);
775 >
776      if(!ffStream->is_open()){
777  
778        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
779 <      ffStream->open( forceFieldFilename.c_str() );
779 >      ffStream->open( forceFieldFilename.c_str(),  
780 >                      ifstream::in | ifstream::binary );
781  
782        //if current directory does not contain the force field file,
783        //try to open it in the path        

Diff Legend

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