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

Comparing trunk/src/UseTheForce/ForceField.cpp (file contents):
Revision 206 by gezelter, Thu Nov 4 20:51:23 2004 UTC vs.
Revision 1195 by cpuglis, Thu Dec 6 20:04:02 2007 UTC

# Line 1 | Line 1
1 + /*
2 + * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 + *
4 + * The University of Notre Dame grants you ("Licensee") a
5 + * non-exclusive, royalty free, license to use, modify and
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
19 + *    notice, this list of conditions and the following disclaimer.
20 + *
21 + * 3. Redistributions in binary form must reproduce the above copyright
22 + *    notice, this list of conditions and the following disclaimer in the
23 + *    documentation and/or other materials provided with the
24 + *    distribution.
25 + *
26 + * This software is provided "AS IS," without a warranty of any
27 + * kind. All express or implied conditions, representations and
28 + * warranties, including any implied warranty of merchantability,
29 + * fitness for a particular purpose or non-infringement, are hereby
30 + * excluded.  The University of Notre Dame and its licensors shall not
31 + * be liable for any damages suffered by licensee as a result of
32 + * using, modifying or distributing the software or its
33 + * derivatives. In no event will the University of Notre Dame or its
34 + * licensors be liable for any lost revenue, profit or data, or for
35 + * direct, indirect, special, consequential, incidental or punitive
36 + * damages, however caused and regardless of the theory of liability,
37 + * arising out of the use of or inability to use software, even if the
38 + * University of Notre Dame has been advised of the possibility of
39 + * such damages.
40 + */
41 +
42 + /**
43 + * @file ForceField.cpp
44 + * @author tlin
45 + * @date 11/04/2004
46 + * @time 22:51am
47 + * @version 1.0
48 + */
49 +  
50   #include "UseTheForce/ForceField.hpp"
51 + #include "utils/simError.h"
52 + #include "UseTheForce/DarkSide/atype_interface.h"
53 + #include "UseTheForce/DarkSide/fForceOptions_interface.h"
54 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
55 + namespace oopse {
56  
57 < AtomType* ForceField::getMatchingAtomType(const string &at) {
57 >  ForceField::ForceField() {
58 >    char* tempPath;
59 >    tempPath = getenv("FORCE_PARAM_PATH");
60  
61 <  map<string, AtomType*>::iterator iter;
62 <  
63 <  iter = atomTypeMap.find(at);
64 <  if (iter != atomTypeMap.end()) {
65 <    return iter->second;
66 <  } else {
11 <    return NULL;
61 >    if (tempPath == NULL) {
62 >      //convert a macro from compiler to a string in c++
63 >      STR_DEFINE(ffPath_, FRC_PATH );
64 >    } else {
65 >      ffPath_ = tempPath;
66 >    }
67    }
13 – }
68  
15 – BondType* ForceField::getMatchingBondType(const string &at1,
16 –                                          const string &at2) {
69  
70 <  map<pair<string,string>, BondType*>::iterator iter;
71 <  vector<BondType*> foundTypes;
70 >  ForceField::~ForceField() {
71 >    deleteAtypes();
72 >    deleteSwitch();
73 >  }
74  
75 <  iter = bondTypeMap.find(pair<at1, at2>);
76 <  if (iter != bondTypeMap.end()) {
77 <    // exact match, so just return it
78 <    return iter->second;
79 <  }
75 >  AtomType* ForceField::getAtomType(const std::string &at) {
76 >    std::vector<std::string> keys;
77 >    keys.push_back(at);
78 >    return atomTypeCont_.find(keys);
79 >  }
80  
81 <  iter = bondTypeMap.find(pair<at2, at1>);
82 <  if (iter != bondTypeMap.end()) {
83 <    // exact match in reverse order, so just return it
84 <    return iter->second;
85 <  }
81 >  BondType* ForceField::getBondType(const std::string &at1,
82 >                                    const std::string &at2) {
83 >    std::vector<std::string> keys;
84 >    keys.push_back(at1);
85 >    keys.push_back(at2);    
86  
87 <  iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>);
88 <  if (iter != bondTypeMap.end()) {
89 <    foundTypes.push_back(iter->second);
87 >    //try exact match first
88 >    BondType* bondType = bondTypeCont_.find(keys);
89 >    if (bondType) {
90 >      return bondType;
91 >    } else {
92 >      //if no exact match found, try wild card match
93 >      return bondTypeCont_.find(keys, wildCardAtomTypeName_);
94 >    }
95    }
96  
97 <  iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>);
98 <  if (iter != bondTypeMap.end()) {
99 <    foundTypes.push_back(iter->second);
97 >  BendType* ForceField::getBendType(const std::string &at1,
98 >                                    const std::string &at2,
99 >                                    const std::string &at3) {
100 >    std::vector<std::string> keys;
101 >    keys.push_back(at1);
102 >    keys.push_back(at2);    
103 >    keys.push_back(at3);    
104 >
105 >    //try exact match first
106 >    BendType* bendType = bendTypeCont_.find(keys);
107 >    if (bendType) {
108 >      return bendType;
109 >    } else {
110 >      //if no exact match found, try wild card match
111 >      return bendTypeCont_.find(keys, wildCardAtomTypeName_);
112 >    }
113    }
114  
115 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
116 <  if (iter != bondTypeMap.end()) {
117 <    foundTypes.push_back(iter->second);
115 >  TorsionType* ForceField::getTorsionType(const std::string &at1,
116 >                                          const std::string &at2,
117 >                                          const std::string &at3,
118 >                                          const std::string &at4) {
119 >    std::vector<std::string> keys;
120 >    keys.push_back(at1);
121 >    keys.push_back(at2);    
122 >    keys.push_back(at3);    
123 >    keys.push_back(at4);    
124 >
125 >    TorsionType* torsionType = torsionTypeCont_.find(keys);
126 >    if (torsionType) {
127 >      return torsionType;
128 >    } else {
129 >      //if no exact match found, try wild card match
130 >      return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
131 >    }
132 >    
133 >    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
134    }
135  
136 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
137 <  if (iter != bondTypeMap.end()) {
138 <    foundTypes.push_back(iter->second);
136 >  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
137 >    std::vector<std::string> keys;
138 >    keys.push_back(at1);
139 >    keys.push_back(at2);    
140 >    
141 >    //try exact match first
142 >    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
143 >    if (nbiType) {
144 >      return nbiType;
145 >    } else {
146 >      //if no exact match found, try wild card match
147 >      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
148 >    }    
149    }
150    
151 <  if (foundTypes.empty()) {
152 <    return NULL;
153 <  } else {
154 <    
151 >  BondType* ForceField::getExactBondType(const std::string &at1,
152 >                                         const std::string &at2){
153 >    std::vector<std::string> keys;
154 >    keys.push_back(at1);
155 >    keys.push_back(at2);    
156 >    return bondTypeCont_.find(keys);
157 >  }
158 >  
159 >  BendType* ForceField::getExactBendType(const std::string &at1,
160 >                                         const std::string &at2,
161 >                                         const std::string &at3){
162 >    std::vector<std::string> keys;
163 >    keys.push_back(at1);
164 >    keys.push_back(at2);    
165 >    keys.push_back(at3);    
166 >    return bendTypeCont_.find(keys);
167 >  }
168 >  
169 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
170 >                                               const std::string &at2,
171 >                                               const std::string &at3,
172 >                                               const std::string &at4){
173 >    std::vector<std::string> keys;
174 >    keys.push_back(at1);
175 >    keys.push_back(at2);    
176 >    keys.push_back(at3);    
177 >    keys.push_back(at4);  
178 >    return torsionTypeCont_.find(keys);
179 >  }
180  
181 <
181 >  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
182 >    std::vector<std::string> keys;
183 >    keys.push_back(at1);
184 >    keys.push_back(at2);    
185 >    return nonBondedInteractionTypeCont_.find(keys);
186 >  }
187  
188  
189 +  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
190 +    std::vector<std::string> keys;
191 +    keys.push_back(at);
192 +    return atomTypeCont_.add(keys, atomType);
193 +  }
194  
195 +  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
196 +                               BondType* bondType) {
197 +    std::vector<std::string> keys;
198 +    keys.push_back(at1);
199 +    keys.push_back(at2);    
200 +    return bondTypeCont_.add(keys, bondType);    
201 +  }
202    
203 +  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
204 +                               const std::string &at3, BendType* bendType) {
205 +    std::vector<std::string> keys;
206 +    keys.push_back(at1);
207 +    keys.push_back(at2);    
208 +    keys.push_back(at3);    
209 +    return bendTypeCont_.add(keys, bendType);
210 +  }
211 +  
212 +  bool ForceField::addTorsionType(const std::string &at1,
213 +                                  const std::string &at2,
214 +                                  const std::string &at3,
215 +                                  const std::string &at4,
216 +                                  TorsionType* torsionType) {
217 +    std::vector<std::string> keys;
218 +    keys.push_back(at1);
219 +    keys.push_back(at2);    
220 +    keys.push_back(at3);    
221 +    keys.push_back(at4);    
222 +    return torsionTypeCont_.add(keys, torsionType);
223 +  }
224  
225 +  bool ForceField::addNonBondedInteractionType(const std::string &at1,
226 +                                               const std::string &at2,
227 +                                               NonBondedInteractionType* nbiType) {
228 +    std::vector<std::string> keys;
229 +    keys.push_back(at1);
230 +    keys.push_back(at2);    
231 +    return nonBondedInteractionTypeCont_.add(keys, nbiType);
232 +  }
233 +  
234 +  RealType ForceField::getRcutFromAtomType(AtomType* at) {
235 +    /**@todo */
236 +    GenericData* data;
237 +    RealType rcut = 0.0;
238 +    
239 +    if (at->isLennardJones()) {
240 +      data = at->getPropertyByName("LennardJones");
241 +      if (data != NULL) {
242 +        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
243 +        
244 +        if (ljData != NULL) {
245 +          LJParam ljParam = ljData->getData();
246 +          
247 +          //by default use 2.5*sigma as cutoff radius
248 +          rcut = 2.5 * ljParam.sigma;
249 +          
250 +        } else {
251 +          sprintf( painCave.errMsg,
252 +                   "Can not cast GenericData to LJParam\n");
253 +          painCave.severity = OOPSE_ERROR;
254 +          painCave.isFatal = 1;
255 +          simError();          
256 +        }            
257 +      } else {
258 +        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
259 +        painCave.severity = OOPSE_ERROR;
260 +        painCave.isFatal = 1;
261 +        simError();          
262 +      }
263 +    }
264 +    return rcut;    
265 +  }
266 +  
267  
268 < BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
269 <                                          const string &at3);
270 < TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
271 <                                                const string &at3, const string &at4);
268 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
269 >    std::string forceFieldFilename(filename);
270 >    ifstrstream* ffStream = new ifstrstream();
271 >    
272 >    //try to open the force filed file in current directory first    
273 >    ffStream->open(forceFieldFilename.c_str());
274 >    if(!ffStream->is_open()){
275  
276 < double ForceField::getRcutForAtomType(AtomType* at);
276 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
277 >      ffStream->open( forceFieldFilename.c_str() );
278  
279 +      //if current directory does not contain the force field file,
280 +      //try to open it in the path        
281 +      if(!ffStream->is_open()){
282  
283 < vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
284 <  
285 <   vector<vector<string> > results;
283 >        sprintf( painCave.errMsg,
284 >                 "Error opening the force field parameter file:\n"
285 >                 "\t%s\n"
286 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
287 >                 "variable?\n",
288 >                 forceFieldFilename.c_str() );
289 >        painCave.severity = OOPSE_ERROR;
290 >        painCave.isFatal = 1;
291 >        simError();
292 >      }
293 >    }  
294 >    return ffStream;
295 >  }
296  
297 <  
298 <
299 <
300 <   vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
301 <    
302 <     int nStrings;
83 <     vector<string> oneResult;
84 <     vector<vector<string> > allResults;
85 <
86 <     nStrings = myAts.size();
87 <
88 <     if (nStrings == 1) {
89 <       oneResult.push_back(wildcardCharacter);
90 <       allResults.push_back(oneResult);
91 <       return allResults;
92 <     } else {
93 <      
94 <       for (i=0; i < nStrings; i++) {
95 <         oneResult = myAts;
96 <         replace(oneResult.begin(), oneResult.end(),
297 >  void ForceField::setFortranForceOptions(){
298 >    ForceOptions theseFortranOptions;
299 >    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
300 >    setfForceOptions(&theseFortranOptions);
301 >  }
302 > } //end namespace oopse

Diff Legend

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