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

Diff Legend

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