ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/heatflux/src/UseTheForce/MnM_FF.cpp
(Generate patch)

Comparing trunk/src/UseTheForce/MnM_FF.cpp (file contents):
Revision 1237 by chuckv, Thu Jul 12 23:21:00 2007 UTC vs.
Revision 1238 by gezelter, Thu Apr 24 21:06:05 2008 UTC

# Line 64 | Line 64
64   #include "UseTheForce/ForceFieldCreator.hpp"
65  
66  
67 < namespace oopse
68 < {
69 <
70 < MnM_FF::MnM_FF()
71 < {
72 <
73 <        //set default force field filename
74 <        setForceFieldFileName("MnM.frc");
75 <
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 < }
106 <
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();
67 > namespace oopse {
68 >  
69 >  MnM_FF::MnM_FF() {    
70 >    
71 >    //set default force field filename
72 >    setForceFieldFileName("MnM.frc");
73 >    
74 >    //The order of adding section parsers is important.
75 >    //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    }
104    
105 <  for (at = atomTypeCont_.beginType(i); at != NULL;
106 <       at = atomTypeCont_.nextType(i)) {
107 <    at->complete();
108 <  }
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 <  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)) {
105 >  void MnM_FF::parse(const std::string& filename) {
106 >    ifstrstream* ffStream;
107 >        
108 >    ffStream = openForceFieldFile(filename);
109      
110 <    /* if we find a metallic atom, we need to compare against
161 <       all other atom types */
110 >    spMan_.parse(*ffStream, *this);
111      
112 <    if (at->isEAM() || at->isSC()) {
112 >    ForceField::AtomTypeContainer::MapTypeIterator i;
113 >    AtomType* at;
114 >    ForceField::AtomTypeContainer::MapTypeIterator j;
115 >    AtomType* at2;
116 >    ForceField::NonBondedInteractionTypeContainer::MapTypeIterator k;
117 >    NonBondedInteractionType* nbit;
118 >    
119 >    for (at = atomTypeCont_.beginType(i); at != NULL;
120 >         at = atomTypeCont_.nextType(i)) {
121 >      at->makeFortranAtomType();
122 >    }
123 >    
124 >    for (at = atomTypeCont_.beginType(i); at != NULL;
125 >         at = atomTypeCont_.nextType(i)) {
126 >      at->complete();
127 >    }
128 >    
129 >    hasSCtypes_ = false;
130 >    for (at = atomTypeCont_.beginType(i); at != NULL;
131 >         at = atomTypeCont_.nextType(i)) {
132 >      if (at->isSC())
133 >        hasSCtypes_ = true;
134 >    }
135 >    
136 >    hasEAMtypes_ = false;
137 >    for (at = atomTypeCont_.beginType(i); at != NULL;
138 >         at = atomTypeCont_.nextType(i)) {
139 >      if (at->isEAM())
140 >        hasEAMtypes_ = true;
141 >    }
142 >    
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 >    /* to handle metal-nonmetal interactions, first we loop over
152 >       all atom types: */
153 >    
154 >    for (at = atomTypeCont_.beginType(i); at != NULL;
155 >         at = atomTypeCont_.nextType(i)) {
156        
157 <      /* loop over all other atom types */
158 <      for (at2 = atomTypeCont_.beginType(j); at2 != NULL;
159 <           at2 = atomTypeCont_.nextType(j)) {
160 <        
161 <        /* if the other partner is not a metallic type, we need to
162 <           look for explicit non-bonded interactions */
163 <        if (!at2->isEAM() && !at2->isSC()) {
157 >      /* if we find a metallic atom, we need to compare against
158 >         all other atom types */
159 >      
160 >      if (at->isEAM() || at->isSC()) {
161 >        
162 >        /* loop over all other atom types */
163 >        for (at2 = atomTypeCont_.beginType(j); at2 != NULL;
164 >             at2 = atomTypeCont_.nextType(j)) {
165 >          
166 >          /* if the other partner is not a metallic type, we need to
167 >             look for explicit non-bonded interactions */
168 >          if (!at2->isEAM() && !at2->isSC()) {
169            
170 <          /* get the name and ident of the metallic atom type */
171 <          std::string at1s = at->getName();
172 <          int atid1 = at->getIdent();
170 >            /* get the name and ident of the metallic atom type */
171 >            std::string at1s = at->getName();
172 >            int atid1 = at->getIdent();
173            
174 <          /* get the name and ident of the nonmetallic atom type */
175 <          std::string at2s = at2->getName();
176 <          int atid2 = at2->getIdent();
174 >            /* get the name and ident of the nonmetallic atom type */
175 >            std::string at2s = at2->getName();
176 >            int atid2 = at2->getIdent();
177            
178 <          /* look for a match in the non-bonded interactions parsed
179 <             from the force field file */
180 <          nbit = getNonBondedInteractionType(at1s, at2s);
178 >            /* look for a match in the non-bonded interactions parsed
179 >               from the force field file */
180 >            nbit = getNonBondedInteractionType(at1s, at2s);
181            
182 <          /* if we found a match (even a partial match), punt to the
183 <             interaction to poke our info down to fortran. */
184 <          if (nbit != NULL)     nbit->tellFortran(atid1, atid2);
185 <        }                      
182 >            /* if we found a match (even a partial match), punt to the
183 >               interaction to poke our info down to fortran. */
184 >            if (nbit != NULL)   nbit->tellFortran(atid1, atid2);
185 >          }                    
186 >        }
187        }
188      }
189 +  
190 +    delete ffStream;
191 +  
192    }
193    
193 –  delete ffStream;
194    
195 < }
196 <
197 <
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) {
206 <
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, "Can not find EAM Parameters\n");
220 <                                painCave.severity = OOPSE_ERROR;
221 <                                painCave.isFatal = 1;
222 <                                simError();
223 <                        }
224 <                }
225 <                else {
226 <                        rcut = ForceField::getRcutFromAtomType(at);
227 <                }
228 <
229 <                return rcut;
230 <        }
231 <
232 <
233 <
234 <
235 <
236 <
237 <        MnM_FF::~MnM_FF() {
238 <                destroyLJTypes();
239 <                destroyStickyTypes();
240 <                if (hasEAMtypes_) destroyEAMTypes();
241 <                if (hasSCtypes_)  destroySCTypes();
242 <        }
195 >  RealType MnM_FF::getRcutFromAtomType(AtomType* at) {
196 >    RealType rcut = 0.0;
197 >    if (at->isEAM()) {
198 >      GenericData* data = at->getPropertyByName("EAM");
199 >      if (data != NULL) {
200 >        EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
201 >        
202 >        if (eamData != NULL) {
203 >          
204 >          EAMParam& eamParam = eamData->getData();
205 >          rcut =  eamParam.rcut;
206 >        }
207 >        else {
208 >          sprintf(painCave.errMsg,
209 >                  "Can not cast GenericData to EAMParam\n");
210 >          painCave.severity = OOPSE_ERROR;
211 >          painCave.isFatal = 1;
212 >          simError();
213 >        }
214 >      }
215 >      else {
216 >        sprintf(painCave.errMsg, "Can not find EAM Parameters\n");
217 >        painCave.severity = OOPSE_ERROR;
218 >        painCave.isFatal = 1;
219 >        simError();
220 >      }
221 >    }
222 >    else {
223 >      rcut = ForceField::getRcutFromAtomType(at);
224 >    }
225 >    
226 >    return rcut;
227 >  }
228 >    
229 >  MnM_FF::~MnM_FF() {
230 >    destroyLJTypes();
231 >    destroyStickyTypes();
232 >    if (hasEAMtypes_) destroyEAMTypes();
233 >    if (hasSCtypes_)  destroySCTypes();
234 >  }
235   } //end namespace oopse
236  

Diff Legend

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