ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/UseTheForce/MnM_FF.cpp
(Generate patch)

Comparing trunk/OOPSE-4/src/UseTheForce/MnM_FF.cpp (file contents):
Revision 3154 by chuckv, Mon Jul 9 21:03:04 2007 UTC vs.
Revision 3164 by chuckv, Thu Jul 12 23:21:00 2007 UTC

# Line 64 | Line 64 | namespace oopse {
64   #include "UseTheForce/ForceFieldCreator.hpp"
65  
66  
67 < namespace oopse {
68 <    
69 <  MnM_FF::MnM_FF(){
67 > namespace oopse
68 > {
69  
70 <    //set default force field filename
71 <    setForceFieldFileName("MnM.frc");
70 > MnM_FF::MnM_FF()
71 > {
72  
73 <    //The order of adding section parsers is important.
74 <    //OptionSectionParser must come first to set options for other parsers
76 <    //DirectionalAtomTypesSectionParser should be added before
77 <    //AtomTypesSectionParser, and these two section parsers will actually
78 <    //create "real" AtomTypes (AtomTypesSectionParser will create AtomType and
79 <    //DirectionalAtomTypesSectionParser will create DirectionalAtomType, which
80 <    //is a subclass of AtomType and should come first). Other AtomTypes Section
81 <    //Parser will not create the "real" AtomType, they only add and set some
82 <    //attribute of the AtomType. Thus their order are not important.
83 <    //AtomTypesSectionParser should be added before other atom type section
84 <    //parsers. Make sure they are added after DirectionalAtomTypesSectionParser
85 <    //and AtomTypesSectionParser. The order of BondTypesSectionParser,
86 <    //BendTypesSectionParser and TorsionTypesSectionParser are not important.
87 <    spMan_.push_back(new OptionSectionParser(forceFieldOptions_));    
88 <    spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
89 <    spMan_.push_back(new AtomTypesSectionParser());
90 <    spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
91 <    spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
92 <    spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
93 <    spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
94 <    spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
95 <    spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
96 <    spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
97 <    spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
98 <                spMan_.push_back(new MetalNonMetalInteractionsSectionParser(forceFieldOptions_));
99 <                spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));    
100 <                spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));    
101 <        spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
102 <    
103 <  }
73 >        //set default force field filename
74 >        setForceFieldFileName("MnM.frc");
75  
76 <  void MnM_FF::parse(const std::string& filename) {
77 <    ifstrstream* ffStream;
78 <        
76 >        //The order of adding section parsers is important.
77 >        //OptionSectionParser must come first to set options for other parsers
78 >        //DirectionalAtomTypesSectionParser should be added before
79 >        //AtomTypesSectionParser, and these two section parsers will actually
80 >        //create "real" AtomTypes (AtomTypesSectionParser will create AtomType and
81 >        //DirectionalAtomTypesSectionParser will create DirectionalAtomType, which
82 >        //is a subclass of AtomType and should come first). Other AtomTypes Section
83 >        //Parser will not create the "real" AtomType, they only add and set some
84 >        //attribute of the AtomType. Thus their order are not important.
85 >        //AtomTypesSectionParser should be added before other atom type section
86 >        //parsers. Make sure they are added after DirectionalAtomTypesSectionParser
87 >        //and AtomTypesSectionParser. The order of BondTypesSectionParser,
88 >        //BendTypesSectionParser and TorsionTypesSectionParser are not important.
89 >        spMan_.push_back(new OptionSectionParser(forceFieldOptions_));
90 >        spMan_.push_back(new DirectionalAtomTypesSectionParser(forceFieldOptions_));
91 >        spMan_.push_back(new AtomTypesSectionParser());
92 >        spMan_.push_back(new LennardJonesAtomTypesSectionParser(forceFieldOptions_));
93 >        spMan_.push_back(new ChargeAtomTypesSectionParser(forceFieldOptions_));
94 >        spMan_.push_back(new MultipoleAtomTypesSectionParser(forceFieldOptions_));
95 >        spMan_.push_back(new StickyAtomTypesSectionParser(forceFieldOptions_));
96 >        spMan_.push_back(new StickyPowerAtomTypesSectionParser(forceFieldOptions_));
97 >        spMan_.push_back(new GayBerneAtomTypesSectionParser(forceFieldOptions_));
98 >        spMan_.push_back(new BondTypesSectionParser(forceFieldOptions_));
99 >        spMan_.push_back(new BendTypesSectionParser(forceFieldOptions_));
100 >        spMan_.push_back(new MetalNonMetalInteractionsSectionParser(forceFieldOptions_));
101 >        spMan_.push_back(new SCAtomTypesSectionParser(forceFieldOptions_));
102 >        spMan_.push_back(new EAMAtomTypesSectionParser(forceFieldOptions_));
103 >        spMan_.push_back(new TorsionTypesSectionParser(forceFieldOptions_));
104  
105 <    ffStream = openForceFieldFile(filename);
105 > }
106  
107 <    spMan_.parse(*ffStream, *this);
107 > void MnM_FF::parse(const std::string& filename) {
108 >  ifstrstream* ffStream;
109 >  
110 >  
111 >  ffStream = openForceFieldFile(filename);
112 >  
113 >  spMan_.parse(*ffStream, *this);
114 >  
115 >  ForceField::AtomTypeContainer::MapTypeIterator i;
116 >  AtomType* at;
117 >  ForceField::AtomTypeContainer::MapTypeIterator j;
118 >  AtomType* at2;
119 >  ForceField::NonBondedInteractionTypeContainer::MapTypeIterator k;
120 >  NonBondedInteractionType* nbit;
121 >  
122 >  for (at = atomTypeCont_.beginType(i); at != NULL;
123 >       at = atomTypeCont_.nextType(i)) {
124 >    at->makeFortranAtomType();
125 >  }
126 >  
127 >  for (at = atomTypeCont_.beginType(i); at != NULL;
128 >       at = atomTypeCont_.nextType(i)) {
129 >    at->complete();
130 >  }
131 >  
132 >  hasSCtypes_ = false;
133 >  for (at = atomTypeCont_.beginType(i); at != NULL;
134 >       at = atomTypeCont_.nextType(i)) {
135 >    if (at->isSC())
136 >      hasSCtypes_ = true;
137 >  }
138  
139 <    ForceField::AtomTypeContainer::MapTypeIterator i;
140 <    AtomType* at;
141 <
142 <    for (at = atomTypeCont_.beginType(i); at != NULL;
143 <         at = atomTypeCont_.nextType(i)) {
144 <      at->makeFortranAtomType();
145 <    }
146 <
147 <    for (at = atomTypeCont_.beginType(i); at != NULL;
148 <         at = atomTypeCont_.nextType(i)) {
149 <      at->complete();
150 <    }
151 <
152 <    hasSCtypes_ = false;
153 <    for (at = atomTypeCont_.beginType(i); at != NULL;
154 <         at = atomTypeCont_.nextType(i)) {
155 <      if (at->isSC())
156 <        hasSCtypes_ = true;
157 <    }
158 <
133 <                hasEAMtypes_ = false;
134 <    for (at = atomTypeCont_.beginType(i); at != NULL;
135 <         at = atomTypeCont_.nextType(i)) {
136 <      if (at->isEAM())
137 <        hasEAMtypes_ = true;
138 <    }    
139 <
139 >  hasEAMtypes_ = false;
140 >  for (at = atomTypeCont_.beginType(i); at != NULL;
141 >       at = atomTypeCont_.nextType(i)) {
142 >    if (at->isEAM())
143 >      hasEAMtypes_ = true;
144 >  }
145 >  
146 >  if (hasEAMtypes_ && hasSCtypes_) {
147 >    sprintf(painCave.errMsg,
148 >            "MnM_FF forcefield cannot use both EAM and Sutton-Chen at the same time\n");
149 >    painCave.severity = OOPSE_ERROR;
150 >    painCave.isFatal = 1;
151 >    simError();
152 >  }
153 >  
154 >  /* to handle metal-nonmetal interactions, first we loop over
155 >     all atom types: */
156 >  
157 >  for (at = atomTypeCont_.beginType(i); at != NULL;
158 >       at = atomTypeCont_.nextType(i)) {
159      
160 <                
161 <
143 <    if (hasEAMtypes_ && hasSCtypes_) {
144 <      sprintf( painCave.errMsg,
145 <                           "MnM_FF forcefield cannot use both EAM and Sutton-Chen at the same time\n");
146 <                  painCave.severity = OOPSE_ERROR;
147 <                  painCave.isFatal = 1;
148 <                  simError();
149 <    }
150 <
151 <    delete ffStream;
160 >    /* if we find a metallic atom, we need to compare against
161 >       all other atom types */
162      
163 +    if (at->isEAM() || at->isSC()) {
164 +      
165 +      /* loop over all other atom types */
166 +      for (at2 = atomTypeCont_.beginType(j); at2 != NULL;
167 +           at2 = atomTypeCont_.nextType(j)) {
168 +        
169 +        /* if the other partner is not a metallic type, we need to
170 +           look for explicit non-bonded interactions */
171 +        if (!at2->isEAM() && !at2->isSC()) {
172 +          
173 +          /* get the name and ident of the metallic atom type */
174 +          std::string at1s = at->getName();
175 +          int atid1 = at->getIdent();
176 +          
177 +          /* get the name and ident of the nonmetallic atom type */
178 +          std::string at2s = at2->getName();
179 +          int atid2 = at2->getIdent();
180 +          
181 +          /* look for a match in the non-bonded interactions parsed
182 +             from the force field file */
183 +          nbit = getNonBondedInteractionType(at1s, at2s);
184 +          
185 +          /* if we found a match (even a partial match), punt to the
186 +             interaction to poke our info down to fortran. */
187 +          if (nbit != NULL)     nbit->tellFortran(atid1, atid2);
188 +        }                      
189 +      }
190 +    }
191    }
192 +  
193 +  delete ffStream;
194 +  
195 + }
196  
197  
198 < RealType MnM_FF::getRcutFromAtomType(AtomType* at)
199 < {
200 <        RealType rcut = 0.0;
201 <        if (at->isEAM()) {
202 <                GenericData* data = at->getPropertyByName("EAM");
203 <                if (data != NULL) {
162 <                        EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
198 >        RealType MnM_FF::getRcutFromAtomType(AtomType* at) {
199 >                RealType rcut = 0.0;
200 >                if (at->isEAM()) {
201 >                        GenericData* data = at->getPropertyByName("EAM");
202 >                        if (data != NULL) {
203 >                                EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
204  
205 <                        if (eamData != NULL) {
205 >                                if (eamData != NULL) {
206  
207 <                                EAMParam& eamParam = eamData->getData();
208 <                                rcut =  eamParam.rcut;
207 >                                        EAMParam& eamParam = eamData->getData();
208 >                                        rcut =  eamParam.rcut;
209 >                                }
210 >                                else {
211 >                                        sprintf(painCave.errMsg,
212 >                                                "Can not cast GenericData to EAMParam\n");
213 >                                        painCave.severity = OOPSE_ERROR;
214 >                                        painCave.isFatal = 1;
215 >                                        simError();
216 >                                }
217                          }
218                          else {
219 <                                sprintf(painCave.errMsg,
171 <                                        "Can not cast GenericData to EAMParam\n");
219 >                                sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
220                                  painCave.severity = OOPSE_ERROR;
221                                  painCave.isFatal = 1;
222                                  simError();
223                          }
224                  }
225                  else {
226 <                        sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
179 <                        painCave.severity = OOPSE_ERROR;
180 <                        painCave.isFatal = 1;
181 <                        simError();
226 >                        rcut = ForceField::getRcutFromAtomType(at);
227                  }
228 +
229 +                return rcut;
230          }
184 –        else {
185 –                rcut = ForceField::getRcutFromAtomType(at);
186 –        }
231  
188 –        return rcut;
189 – }
232  
233  
234  
235  
236  
237 <
238 <  MnM_FF::~MnM_FF(){
239 <    destroyLJTypes();
198 <    destroyStickyTypes();
237 >        MnM_FF::~MnM_FF() {
238 >                destroyLJTypes();
239 >                destroyStickyTypes();
240                  if (hasEAMtypes_) destroyEAMTypes();
241                  if (hasSCtypes_)  destroySCTypes();
242 <  }
242 >        }
243   } //end namespace oopse
244 +

Diff Legend

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