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 939 by gezelter, Thu Apr 20 18:24:24 2006 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;
31 <  }
81 >  BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
82 >    std::vector<std::string> keys;
83 >    keys.push_back(at1);
84 >    keys.push_back(at2);    
85  
86 <  iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>);
87 <  if (iter != bondTypeMap.end()) {
88 <    foundTypes.push_back(iter->second);
86 >    //try exact match first
87 >    BondType* bondType = bondTypeCont_.find(keys);
88 >    if (bondType) {
89 >      return bondType;
90 >    } else {
91 >      //if no exact match found, try wild card match
92 >      return bondTypeCont_.find(keys, wildCardAtomTypeName_);
93 >    }
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, const std::string &at2,
98 >                                    const std::string &at3) {
99 >    std::vector<std::string> keys;
100 >    keys.push_back(at1);
101 >    keys.push_back(at2);    
102 >    keys.push_back(at3);    
103 >
104 >    //try exact match first
105 >    BendType* bendType = bendTypeCont_.find(keys);
106 >    if (bendType) {
107 >      return bendType;
108 >    } else {
109 >      //if no exact match found, try wild card match
110 >      return bendTypeCont_.find(keys, wildCardAtomTypeName_);
111 >    }
112    }
113  
114 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
115 <  if (iter != bondTypeMap.end()) {
116 <    foundTypes.push_back(iter->second);
114 >  TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
115 >                                          const std::string &at3, const std::string &at4) {
116 >    std::vector<std::string> keys;
117 >    keys.push_back(at1);
118 >    keys.push_back(at2);    
119 >    keys.push_back(at3);    
120 >    keys.push_back(at4);    
121 >
122 >    TorsionType* torsionType = torsionTypeCont_.find(keys);
123 >    if (torsionType) {
124 >      return torsionType;
125 >    } else {
126 >      //if no exact match found, try wild card match
127 >      return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
128 >    }
129 >    
130 >    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
131 >
132    }
133  
134 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
135 <  if (iter != bondTypeMap.end()) {
136 <    foundTypes.push_back(iter->second);
134 >  BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
135 >    std::vector<std::string> keys;
136 >    keys.push_back(at1);
137 >    keys.push_back(at2);    
138 >    return bondTypeCont_.find(keys);
139    }
52 –  
53 –  if (foundTypes.empty()) {
54 –    return NULL;
55 –  } else {
56 –    
140  
141 <
141 >  BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
142 >                                         const std::string &at3){
143 >    std::vector<std::string> keys;
144 >    keys.push_back(at1);
145 >    keys.push_back(at2);    
146 >    keys.push_back(at3);    
147 >    return bendTypeCont_.find(keys);
148 >  }
149  
150 +  TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
151 +                                               const std::string &at3, const std::string &at4){
152 +    std::vector<std::string> keys;
153 +    keys.push_back(at1);
154 +    keys.push_back(at2);    
155 +    keys.push_back(at3);    
156 +    keys.push_back(at4);  
157 +    return torsionTypeCont_.find(keys);
158 +  }
159 +  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
160 +    std::vector<std::string> keys;
161 +    keys.push_back(at);
162 +    return atomTypeCont_.add(keys, atomType);
163 +  }
164  
165 +  bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
166 +    std::vector<std::string> keys;
167 +    keys.push_back(at1);
168 +    keys.push_back(at2);    
169 +    return bondTypeCont_.add(keys, bondType);
170  
171 <  
171 >  }
172  
173 +  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
174 +                               const std::string &at3, BendType* bendType) {
175 +    std::vector<std::string> keys;
176 +    keys.push_back(at1);
177 +    keys.push_back(at2);    
178 +    keys.push_back(at3);    
179 +    return bendTypeCont_.add(keys, bendType);
180 +  }
181  
182 < BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
183 <                                          const string &at3);
184 < TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
185 <                                                const string &at3, const string &at4);
182 >  bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
183 >                                  const std::string &at3, const std::string &at4, TorsionType* torsionType) {
184 >    std::vector<std::string> keys;
185 >    keys.push_back(at1);
186 >    keys.push_back(at2);    
187 >    keys.push_back(at3);    
188 >    keys.push_back(at4);    
189 >    return torsionTypeCont_.add(keys, torsionType);
190 >  }
191  
192 < double ForceField::getRcutForAtomType(AtomType* at);
192 >  double ForceField::getRcutFromAtomType(AtomType* at) {
193 >    /**@todo */
194 >    GenericData* data;
195 >    double rcut = 0.0;
196  
197 +    if (at->isLennardJones()) {
198 +      data = at->getPropertyByName("LennardJones");
199 +      if (data != NULL) {
200 +        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
201  
202 < vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
203 <  
75 <   vector<vector<string> > results;
202 >        if (ljData != NULL) {
203 >          LJParam ljParam = ljData->getData();
204  
205 <  
205 >          //by default use 2.5*sigma as cutoff radius
206 >          rcut = 2.5 * ljParam.sigma;
207 >                
208 >        } else {
209 >          sprintf( painCave.errMsg,
210 >                   "Can not cast GenericData to LJParam\n");
211 >          painCave.severity = OOPSE_ERROR;
212 >          painCave.isFatal = 1;
213 >          simError();          
214 >        }            
215 >      } else {
216 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
217 >        painCave.severity = OOPSE_ERROR;
218 >        painCave.isFatal = 1;
219 >        simError();          
220 >      }
221 >    }
222  
223 +    return rcut;    
224 +  }
225  
80 –   vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
81 –    
82 –     int nStrings;
83 –     vector<string> oneResult;
84 –     vector<vector<string> > allResults;
226  
227 <     nStrings = myAts.size();
227 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
228 >    std::string forceFieldFilename(filename);
229 >    ifstrstream* ffStream = new ifstrstream();
230 >    
231 >    //try to open the force filed file in current directory first    
232 >    ffStream->open(forceFieldFilename.c_str());
233 >    if(!ffStream->is_open()){
234  
235 <     if (nStrings == 1) {
236 <       oneResult.push_back(wildcardCharacter);
237 <       allResults.push_back(oneResult);
238 <       return allResults;
239 <     } else {
240 <      
241 <       for (i=0; i < nStrings; i++) {
242 <         oneResult = myAts;
243 <         replace(oneResult.begin(), oneResult.end(),
235 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
236 >      ffStream->open( forceFieldFilename.c_str() );
237 >
238 >      //if current directory does not contain the force field file,
239 >      //try to open it in the path        
240 >      if(!ffStream->is_open()){
241 >
242 >        sprintf( painCave.errMsg,
243 >                 "Error opening the force field parameter file:\n"
244 >                 "\t%s\n"
245 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
246 >                 "variable?\n",
247 >                 forceFieldFilename.c_str() );
248 >        painCave.severity = OOPSE_ERROR;
249 >        painCave.isFatal = 1;
250 >        simError();
251 >      }
252 >    }  
253 >
254 >    return ffStream;
255 >
256 >  }
257 >
258 >  void ForceField::setFortranForceOptions(){
259 >    ForceOptions theseFortranOptions;
260 >    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
261 >    setfForceOptions(&theseFortranOptions);
262 >  }
263 > } //end namespace oopse

Diff Legend

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