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 1282 by gezelter, Wed Jul 30 18:11:19 2008 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 <algorithm>
51   #include "UseTheForce/ForceField.hpp"
52 + #include "utils/simError.h"
53 + #include "utils/Tuple.hpp"
54 + #include "UseTheForce/DarkSide/atype_interface.h"
55 + #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
57 + namespace oopse {
58  
59 < AtomType* ForceField::getMatchingAtomType(const string &at) {
59 >  ForceField::ForceField() {
60 >    char* tempPath;
61 >    tempPath = getenv("FORCE_PARAM_PATH");
62  
63 <  map<string, AtomType*>::iterator iter;
64 <  
65 <  iter = atomTypeMap.find(at);
66 <  if (iter != atomTypeMap.end()) {
67 <    return iter->second;
68 <  } else {
11 <    return NULL;
63 >    if (tempPath == NULL) {
64 >      //convert a macro from compiler to a string in c++
65 >      STR_DEFINE(ffPath_, FRC_PATH );
66 >    } else {
67 >      ffPath_ = tempPath;
68 >    }
69    }
13 – }
70  
15 – BondType* ForceField::getMatchingBondType(const string &at1,
16 –                                          const string &at2) {
71  
72 <  map<pair<string,string>, BondType*>::iterator iter;
73 <  vector<BondType*> foundTypes;
72 >  ForceField::~ForceField() {
73 >    deleteAtypes();
74 >    deleteSwitch();
75 >  }
76  
77 <  iter = bondTypeMap.find(pair<at1, at2>);
78 <  if (iter != bondTypeMap.end()) {
79 <    // exact match, so just return it
80 <    return iter->second;
81 <  }
77 >  AtomType* ForceField::getAtomType(const std::string &at) {
78 >    std::vector<std::string> keys;
79 >    keys.push_back(at);
80 >    return atomTypeCont_.find(keys);
81 >  }
82  
83 <  iter = bondTypeMap.find(pair<at2, at1>);
84 <  if (iter != bondTypeMap.end()) {
85 <    // exact match in reverse order, so just return it
86 <    return iter->second;
87 <  }
83 >  BondType* ForceField::getBondType(const std::string &at1,
84 >                                    const std::string &at2) {
85 >    std::vector<std::string> keys;
86 >    keys.push_back(at1);
87 >    keys.push_back(at2);    
88  
89 <  iter = bondTypeMap.find(pair<at1, wildCardAtomTypeName>);
90 <  if (iter != bondTypeMap.end()) {
91 <    foundTypes.push_back(iter->second);
92 <  }
89 >    //try exact match first
90 >    BondType* bondType = bondTypeCont_.find(keys);
91 >    if (bondType) {
92 >      return bondType;
93 >    } else {
94 >      AtomType* atype1;
95 >      AtomType* atype2;
96 >      std::vector<std::string> at1key;
97 >      at1key.push_back(at1);
98 >      atype1 = atomTypeCont_.find(at1key);
99 >  
100 >      std::vector<std::string> at2key;
101 >      at2key.push_back(at2);
102 >      atype2 = atomTypeCont_.find(at2key);
103  
104 <  iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>);
105 <  if (iter != bondTypeMap.end()) {
106 <    foundTypes.push_back(iter->second);
41 <  }
104 >      // query atom types for their chains of responsibility
105 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
106 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
107  
108 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
109 <  if (iter != bondTypeMap.end()) {
45 <    foundTypes.push_back(iter->second);
46 <  }
108 >      std::vector<AtomType*>::iterator i;
109 >      std::vector<AtomType*>::iterator j;
110  
111 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
112 <  if (iter != bondTypeMap.end()) {
113 <    foundTypes.push_back(iter->second);
111 >      int ii = 0;
112 >      int jj = 0;
113 >      int bondTypeScore;
114 >
115 >      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
116 >
117 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
118 >        jj = 0;
119 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
120 >
121 >          bondTypeScore = ii + jj;
122 >
123 >          std::vector<std::string> myKeys;
124 >          myKeys.push_back((*i)->getName());
125 >          myKeys.push_back((*j)->getName());
126 >
127 >          std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
128 >          BondType* bondType = bondTypeCont_.find(myKeys);
129 >          if (bondType) {
130 >            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
131 >          }
132 >          jj++;
133 >        }
134 >        ii++;
135 >      }
136 >
137 >
138 >      if (foundBonds.size() > 0) {
139 >        // sort the foundBonds by the score:
140 >        std::sort(foundBonds.begin(), foundBonds.end());
141 >    
142 >        int bestScore = foundBonds[0].first;
143 >        std::vector<std::string> theKeys = foundBonds[0].second;
144 >        
145 >        std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
146 >        BondType* bestType = bondTypeCont_.find(theKeys);
147 >        
148 >        return bestType;
149 >      } else {
150 >        //if no exact match found, try wild card match
151 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
152 >      }
153 >    }
154    }
155    
156 <  if (foundTypes.empty()) {
157 <    return NULL;
158 <  } else {
159 <    
156 >  BendType* ForceField::getBendType(const std::string &at1,
157 >                                    const std::string &at2,
158 >                                    const std::string &at3) {
159 >    std::vector<std::string> keys;
160 >    keys.push_back(at1);
161 >    keys.push_back(at2);    
162 >    keys.push_back(at3);    
163  
164 <
164 >    //try exact match first
165 >    BendType* bendType = bendTypeCont_.find(keys);
166 >    if (bendType) {
167 >      return bendType;
168 >    } else {
169  
170 +      AtomType* atype1;
171 +      AtomType* atype2;
172 +      AtomType* atype3;
173 +      std::vector<std::string> at1key;
174 +      at1key.push_back(at1);
175 +      atype1 = atomTypeCont_.find(at1key);
176 +  
177 +      std::vector<std::string> at2key;
178 +      at2key.push_back(at2);
179 +      atype2 = atomTypeCont_.find(at2key);
180  
181 +      std::vector<std::string> at3key;
182 +      at3key.push_back(at3);
183 +      atype3 = atomTypeCont_.find(at3key);
184  
185 <  
185 >      // query atom types for their chains of responsibility
186 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
187 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
188 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
189  
190 +      std::vector<AtomType*>::iterator i;
191 +      std::vector<AtomType*>::iterator j;
192 +      std::vector<AtomType*>::iterator k;
193  
194 < BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
195 <                                          const string &at3);
196 < TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
197 <                                                const string &at3, const string &at4);
194 >      int ii = 0;
195 >      int jj = 0;
196 >      int kk = 0;
197 >      int IKscore;
198  
199 < double ForceField::getRcutForAtomType(AtomType* at);
199 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
200  
201 +      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
202 +        ii = 0;
203 +        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
204 +          kk = 0;
205 +          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
206 +          
207 +            IKscore = ii + kk;
208  
209 < vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
210 <  
211 <   vector<vector<string> > results;
209 >            std::vector<std::string> myKeys;
210 >            myKeys.push_back((*i)->getName());
211 >            myKeys.push_back((*j)->getName());
212 >            myKeys.push_back((*k)->getName());
213  
214 <  
214 >            BendType* bendType = bendTypeCont_.find(myKeys);
215 >            if (bendType) {
216 >              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
217 >            }
218 >            kk++;
219 >          }
220 >          ii++;
221 >        }
222 >        jj++;
223 >      }
224 >      
225 >      if (foundBends.size() > 0) {
226 >        std::sort(foundBends.begin(), foundBends.end());
227 >        int jscore = foundBends[0].first;
228 >        int ikscore = foundBends[0].second;
229 >        std::vector<std::string> theKeys = foundBends[0].third;
230 >        
231 >        std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
232 >        
233 >        BendType* bestType = bendTypeCont_.find(theKeys);  
234 >        return bestType;
235 >      } else {        
236 >        //if no exact match found, try wild card match
237 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
238 >      }
239 >    }
240 >  }
241  
242 +  TorsionType* ForceField::getTorsionType(const std::string &at1,
243 +                                          const std::string &at2,
244 +                                          const std::string &at3,
245 +                                          const std::string &at4) {
246 +    std::vector<std::string> keys;
247 +    keys.push_back(at1);
248 +    keys.push_back(at2);    
249 +    keys.push_back(at3);    
250 +    keys.push_back(at4);    
251  
80 –   vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
81 –    
82 –     int nStrings;
83 –     vector<string> oneResult;
84 –     vector<vector<string> > allResults;
252  
253 <     nStrings = myAts.size();
253 >    //try exact match first
254 >    TorsionType* torsionType = torsionTypeCont_.find(keys);
255 >    if (torsionType) {
256 >      return torsionType;
257 >    } else {
258  
259 <     if (nStrings == 1) {
260 <       oneResult.push_back(wildcardCharacter);
261 <       allResults.push_back(oneResult);
262 <       return allResults;
263 <     } else {
264 <      
265 <       for (i=0; i < nStrings; i++) {
266 <         oneResult = myAts;
267 <         replace(oneResult.begin(), oneResult.end(),
259 >      AtomType* atype1;
260 >      AtomType* atype2;
261 >      AtomType* atype3;
262 >      AtomType* atype4;
263 >      std::vector<std::string> at1key;
264 >      at1key.push_back(at1);
265 >      atype1 = atomTypeCont_.find(at1key);
266 >  
267 >      std::vector<std::string> at2key;
268 >      at2key.push_back(at2);
269 >      atype2 = atomTypeCont_.find(at2key);
270 >
271 >      std::vector<std::string> at3key;
272 >      at3key.push_back(at3);
273 >      atype3 = atomTypeCont_.find(at3key);
274 >
275 >      std::vector<std::string> at4key;
276 >      at4key.push_back(at4);
277 >      atype4 = atomTypeCont_.find(at4key);
278 >
279 >      // query atom types for their chains of responsibility
280 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
281 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
282 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
283 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
284 >
285 >      std::vector<AtomType*>::iterator i;
286 >      std::vector<AtomType*>::iterator j;
287 >      std::vector<AtomType*>::iterator k;
288 >      std::vector<AtomType*>::iterator l;
289 >
290 >      int ii = 0;
291 >      int jj = 0;
292 >      int kk = 0;
293 >      int ll = 0;
294 >      int ILscore;
295 >      int JKscore;
296 >
297 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
298 >
299 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
300 >        kk = 0;
301 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
302 >          ii = 0;      
303 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
304 >            ll = 0;
305 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
306 >          
307 >              ILscore = ii + ll;
308 >              JKscore = jj + kk;
309 >
310 >              std::vector<std::string> myKeys;
311 >              myKeys.push_back((*i)->getName());
312 >              myKeys.push_back((*j)->getName());
313 >              myKeys.push_back((*k)->getName());
314 >              myKeys.push_back((*l)->getName());
315 >
316 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
317 >              if (torsionType) {
318 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
319 >              }
320 >              ll++;
321 >            }
322 >            ii++;
323 >          }
324 >          kk++;
325 >        }
326 >        jj++;
327 >      }
328 >      
329 >      if (foundTorsions.size() > 0) {
330 >        std::sort(foundTorsions.begin(), foundTorsions.end());
331 >        int jkscore = foundTorsions[0].first;
332 >        int ilscore = foundTorsions[0].second;
333 >        std::vector<std::string> theKeys = foundTorsions[0].third;
334 >        
335 >        std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
336 >                
337 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
338 >        return bestType;
339 >      } else {
340 >        //if no exact match found, try wild card match
341 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
342 >      }
343 >    }
344 >  }
345 >
346 >  InversionType* ForceField::getInversionType(const std::string &at1,
347 >                                              const std::string &at2,
348 >                                              const std::string &at3,
349 >                                              const std::string &at4) {
350 >    std::vector<std::string> keys;
351 >    keys.push_back(at1);
352 >    keys.push_back(at2);    
353 >    keys.push_back(at3);    
354 >    keys.push_back(at4);    
355 >
356 >    //try exact match first
357 >    InversionType* inversionType = inversionTypeCont_.find(keys);
358 >    if (inversionType) {
359 >      return inversionType;
360 >    } else {
361 >      
362 >      AtomType* atype1;
363 >      AtomType* atype2;
364 >      AtomType* atype3;
365 >      AtomType* atype4;
366 >      std::vector<std::string> at1key;
367 >      at1key.push_back(at1);
368 >      atype1 = atomTypeCont_.find(at1key);
369 >      
370 >      std::vector<std::string> at2key;
371 >      at2key.push_back(at2);
372 >      atype2 = atomTypeCont_.find(at2key);
373 >      
374 >      std::vector<std::string> at3key;
375 >      at3key.push_back(at3);
376 >      atype3 = atomTypeCont_.find(at3key);
377 >      
378 >      std::vector<std::string> at4key;
379 >      at4key.push_back(at4);
380 >      atype4 = atomTypeCont_.find(at4key);
381 >
382 >      // query atom types for their chains of responsibility
383 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
384 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
385 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
386 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
387 >
388 >      std::vector<AtomType*>::iterator i;
389 >      std::vector<AtomType*>::iterator j;
390 >      std::vector<AtomType*>::iterator k;
391 >      std::vector<AtomType*>::iterator l;
392 >
393 >      int ii = 0;
394 >      int jj = 0;
395 >      int kk = 0;
396 >      int ll = 0;
397 >      int Iscore;
398 >      int JKLscore;
399 >      
400 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
401 >      
402 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
403 >        kk = 0;
404 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
405 >          ii = 0;      
406 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
407 >            ll = 0;
408 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
409 >              
410 >              Iscore = ii;
411 >              JKLscore = jj + kk + ll;
412 >              
413 >              std::vector<std::string> myKeys;
414 >              myKeys.push_back((*i)->getName());
415 >              myKeys.push_back((*j)->getName());
416 >              myKeys.push_back((*k)->getName());
417 >              myKeys.push_back((*l)->getName());
418 >              
419 >              InversionType* inversionType = inversionTypeCont_.find(myKeys);
420 >              if (inversionType) {
421 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
422 >              }
423 >              ll++;
424 >            }
425 >            ii++;
426 >          }
427 >          kk++;
428 >        }
429 >        jj++;
430 >      }
431 >        
432 >      if (foundInversions.size() > 0) {
433 >        std::sort(foundInversions.begin(), foundInversions.end());
434 >        int iscore = foundInversions[0].first;
435 >        int jklscore = foundInversions[0].second;
436 >        std::vector<std::string> theKeys = foundInversions[0].third;
437 >        
438 >        std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
439 >                
440 >        InversionType* bestType = inversionTypeCont_.find(theKeys);
441 >        return bestType;
442 >      } else {
443 >        //if no exact match found, try wild card match
444 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
445 >      }
446 >    }
447 >  }
448 >  
449 >  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
450 >    std::vector<std::string> keys;
451 >    keys.push_back(at1);
452 >    keys.push_back(at2);    
453 >    
454 >    //try exact match first
455 >    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
456 >    if (nbiType) {
457 >      return nbiType;
458 >    } else {
459 >      //if no exact match found, try wild card match
460 >      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
461 >    }    
462 >  }
463 >  
464 >  BondType* ForceField::getExactBondType(const std::string &at1,
465 >                                         const std::string &at2){
466 >    std::vector<std::string> keys;
467 >    keys.push_back(at1);
468 >    keys.push_back(at2);    
469 >    return bondTypeCont_.find(keys);
470 >  }
471 >  
472 >  BendType* ForceField::getExactBendType(const std::string &at1,
473 >                                         const std::string &at2,
474 >                                         const std::string &at3){
475 >    std::vector<std::string> keys;
476 >    keys.push_back(at1);
477 >    keys.push_back(at2);    
478 >    keys.push_back(at3);    
479 >    return bendTypeCont_.find(keys);
480 >  }
481 >  
482 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
483 >                                               const std::string &at2,
484 >                                               const std::string &at3,
485 >                                               const std::string &at4){
486 >    std::vector<std::string> keys;
487 >    keys.push_back(at1);
488 >    keys.push_back(at2);    
489 >    keys.push_back(at3);    
490 >    keys.push_back(at4);  
491 >    return torsionTypeCont_.find(keys);
492 >  }
493 >  
494 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
495 >                                                   const std::string &at2,
496 >                                                   const std::string &at3,
497 >                                                   const std::string &at4){
498 >    std::vector<std::string> keys;
499 >    keys.push_back(at1);
500 >    keys.push_back(at2);    
501 >    keys.push_back(at3);    
502 >    keys.push_back(at4);  
503 >    return inversionTypeCont_.find(keys);
504 >  }
505 >  
506 >  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
507 >    std::vector<std::string> keys;
508 >    keys.push_back(at1);
509 >    keys.push_back(at2);    
510 >    return nonBondedInteractionTypeCont_.find(keys);
511 >  }
512 >  
513 >
514 >  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
515 >    std::vector<std::string> keys;
516 >    keys.push_back(at);
517 >    return atomTypeCont_.add(keys, atomType);
518 >  }
519 >
520 >  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
521 >    std::vector<std::string> keys;
522 >    keys.push_back(at);
523 >    return atomTypeCont_.replace(keys, atomType);
524 >  }
525 >
526 >  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
527 >                               BondType* bondType) {
528 >    std::vector<std::string> keys;
529 >    keys.push_back(at1);
530 >    keys.push_back(at2);    
531 >    return bondTypeCont_.add(keys, bondType);    
532 >  }
533 >  
534 >  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
535 >                               const std::string &at3, BendType* bendType) {
536 >    std::vector<std::string> keys;
537 >    keys.push_back(at1);
538 >    keys.push_back(at2);    
539 >    keys.push_back(at3);    
540 >    return bendTypeCont_.add(keys, bendType);
541 >  }
542 >  
543 >  bool ForceField::addTorsionType(const std::string &at1,
544 >                                  const std::string &at2,
545 >                                  const std::string &at3,
546 >                                  const std::string &at4,
547 >                                  TorsionType* torsionType) {
548 >    std::vector<std::string> keys;
549 >    keys.push_back(at1);
550 >    keys.push_back(at2);    
551 >    keys.push_back(at3);    
552 >    keys.push_back(at4);    
553 >    return torsionTypeCont_.add(keys, torsionType);
554 >  }
555 >
556 >  bool ForceField::addInversionType(const std::string &at1,
557 >                                    const std::string &at2,
558 >                                    const std::string &at3,
559 >                                    const std::string &at4,
560 >                                    InversionType* inversionType) {
561 >    std::vector<std::string> keys;
562 >    keys.push_back(at1);
563 >    keys.push_back(at2);    
564 >    keys.push_back(at3);    
565 >    keys.push_back(at4);    
566 >    return inversionTypeCont_.add(keys, inversionType);
567 >  }
568 >  
569 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
570 >                                               const std::string &at2,
571 >                                               NonBondedInteractionType* nbiType) {
572 >    std::vector<std::string> keys;
573 >    keys.push_back(at1);
574 >    keys.push_back(at2);    
575 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
576 >  }
577 >  
578 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
579 >    /**@todo */
580 >    GenericData* data;
581 >    RealType rcut = 0.0;
582 >    
583 >    if (at->isLennardJones()) {
584 >      data = at->getPropertyByName("LennardJones");
585 >      if (data != NULL) {
586 >        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
587 >        
588 >        if (ljData != NULL) {
589 >          LJParam ljParam = ljData->getData();
590 >          
591 >          //by default use 2.5*sigma as cutoff radius
592 >          rcut = 2.5 * ljParam.sigma;
593 >          
594 >        } else {
595 >          sprintf( painCave.errMsg,
596 >                   "Can not cast GenericData to LJParam\n");
597 >          painCave.severity = OOPSE_ERROR;
598 >          painCave.isFatal = 1;
599 >          simError();          
600 >        }            
601 >      } else {
602 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
603 >        painCave.severity = OOPSE_ERROR;
604 >        painCave.isFatal = 1;
605 >        simError();          
606 >      }
607 >    }
608 >    return rcut;    
609 >  }
610 >  
611 >
612 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
613 >    std::string forceFieldFilename(filename);
614 >    ifstrstream* ffStream = new ifstrstream();
615 >    
616 >    //try to open the force filed file in current directory first    
617 >    ffStream->open(forceFieldFilename.c_str());
618 >    if(!ffStream->is_open()){
619 >
620 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
621 >      ffStream->open( forceFieldFilename.c_str() );
622 >
623 >      //if current directory does not contain the force field file,
624 >      //try to open it in the path        
625 >      if(!ffStream->is_open()){
626 >
627 >        sprintf( painCave.errMsg,
628 >                 "Error opening the force field parameter file:\n"
629 >                 "\t%s\n"
630 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
631 >                 "variable?\n",
632 >                 forceFieldFilename.c_str() );
633 >        painCave.severity = OOPSE_ERROR;
634 >        painCave.isFatal = 1;
635 >        simError();
636 >      }
637 >    }  
638 >    return ffStream;
639 >  }
640 >
641 >  void ForceField::setFortranForceOptions(){
642 >    ForceOptions theseFortranOptions;
643 >    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
644 >    setfForceOptions(&theseFortranOptions);
645 >  }
646 > } //end namespace oopse

Diff Legend

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