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.
branches/development/src/UseTheForce/ForceField.cpp (file contents), Revision 1530 by gezelter, Tue Dec 28 21:47:55 2010 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. Redistributions of source code must retain the above copyright
10 + *    notice, this list of conditions and the following disclaimer.
11 + *
12 + * 2. Redistributions in binary form must reproduce the above copyright
13 + *    notice, this list of conditions and the following disclaimer in the
14 + *    documentation and/or other materials provided with the
15 + *    distribution.
16 + *
17 + * This software is provided "AS IS," without a warranty of any
18 + * kind. All express or implied conditions, representations and
19 + * warranties, including any implied warranty of merchantability,
20 + * fitness for a particular purpose or non-infringement, are hereby
21 + * excluded.  The University of Notre Dame and its licensors shall not
22 + * be liable for any damages suffered by licensee as a result of
23 + * using, modifying or distributing the software or its
24 + * derivatives. In no event will the University of Notre Dame or its
25 + * licensors be liable for any lost revenue, profit or data, or for
26 + * direct, indirect, special, consequential, incidental or punitive
27 + * damages, however caused and regardless of the theory of liability,
28 + * arising out of the use of or inability to use software, even if the
29 + * University of Notre Dame has been advised of the possibility of
30 + * such damages.
31 + *
32 + * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33 + * research, please cite the appropriate papers when you publish your
34 + * work.  Good starting points are:
35 + *                                                                      
36 + * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37 + * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 + * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 + * [4]  Vardeman & Gezelter, in progress (2009).                        
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 + namespace OpenMD {
56  
57 < AtomType* ForceField::getMatchingAtomType(const string &at) {
57 >  ForceField::ForceField() {
58  
59 <  map<string, AtomType*>::iterator iter;
60 <  
61 <  iter = atomTypeMap.find(at);
62 <  if (iter != atomTypeMap.end()) {
63 <    return iter->second;
64 <  } else {
65 <    return NULL;
59 >    char* tempPath;
60 >    tempPath = getenv("FORCE_PARAM_PATH");
61 >    
62 >    if (tempPath == NULL) {
63 >      //convert a macro from compiler to a string in c++
64 >      STR_DEFINE(ffPath_, FRC_PATH );
65 >    } else {
66 >      ffPath_ = tempPath;
67 >    }
68    }
13 }
69  
15 BondType* ForceField::getMatchingBondType(const string &at1,
16                                          const string &at2) {
70  
71 <  map<pair<string,string>, BondType*>::iterator iter;
72 <  vector<BondType*> foundTypes;
71 >  ForceField::~ForceField() {
72 >    deleteAtypes();
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);
90 <  }
87 >    //try exact match first
88 >    BondType* bondType = bondTypeCont_.find(keys);
89 >    if (bondType) {
90 >      return bondType;
91 >    } else {
92 >      AtomType* atype1;
93 >      AtomType* atype2;
94 >      std::vector<std::string> at1key;
95 >      at1key.push_back(at1);
96 >      atype1 = atomTypeCont_.find(at1key);
97 >  
98 >      std::vector<std::string> at2key;
99 >      at2key.push_back(at2);
100 >      atype2 = atomTypeCont_.find(at2key);
101  
102 <  iter = bondTypeMap.find(pair<at2, wildCardAtomTypeName>);
103 <  if (iter != bondTypeMap.end()) {
104 <    foundTypes.push_back(iter->second);
41 <  }
102 >      // query atom types for their chains of responsibility
103 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
104 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
105  
106 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at1>);
107 <  if (iter != bondTypeMap.end()) {
45 <    foundTypes.push_back(iter->second);
46 <  }
106 >      std::vector<AtomType*>::iterator i;
107 >      std::vector<AtomType*>::iterator j;
108  
109 <  iter = bondTypeMap.find(pair<wildCardAtomTypeName, at2>);
110 <  if (iter != bondTypeMap.end()) {
111 <    foundTypes.push_back(iter->second);
109 >      int ii = 0;
110 >      int jj = 0;
111 >      int bondTypeScore;
112 >
113 >      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
114 >
115 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
116 >        jj = 0;
117 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
118 >
119 >          bondTypeScore = ii + jj;
120 >
121 >          std::vector<std::string> myKeys;
122 >          myKeys.push_back((*i)->getName());
123 >          myKeys.push_back((*j)->getName());
124 >
125 >          BondType* bondType = bondTypeCont_.find(myKeys);
126 >          if (bondType) {
127 >            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
128 >          }
129 >          jj++;
130 >        }
131 >        ii++;
132 >      }
133 >
134 >
135 >      if (foundBonds.size() > 0) {
136 >        // sort the foundBonds by the score:
137 >        std::sort(foundBonds.begin(), foundBonds.end());
138 >    
139 >        int bestScore = foundBonds[0].first;
140 >        std::vector<std::string> theKeys = foundBonds[0].second;
141 >        
142 >        BondType* bestType = bondTypeCont_.find(theKeys);
143 >        
144 >        return bestType;
145 >      } else {
146 >        //if no exact match found, try wild card match
147 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
148 >      }
149 >    }
150    }
151    
152 <  if (foundTypes.empty()) {
153 <    return NULL;
154 <  } else {
155 <    
152 >  BendType* ForceField::getBendType(const std::string &at1,
153 >                                    const std::string &at2,
154 >                                    const std::string &at3) {
155 >    std::vector<std::string> keys;
156 >    keys.push_back(at1);
157 >    keys.push_back(at2);    
158 >    keys.push_back(at3);    
159  
160 <
160 >    //try exact match first
161 >    BendType* bendType = bendTypeCont_.find(keys);
162 >    if (bendType) {
163 >      return bendType;
164 >    } else {
165  
166 +      AtomType* atype1;
167 +      AtomType* atype2;
168 +      AtomType* atype3;
169 +      std::vector<std::string> at1key;
170 +      at1key.push_back(at1);
171 +      atype1 = atomTypeCont_.find(at1key);
172 +  
173 +      std::vector<std::string> at2key;
174 +      at2key.push_back(at2);
175 +      atype2 = atomTypeCont_.find(at2key);
176  
177 +      std::vector<std::string> at3key;
178 +      at3key.push_back(at3);
179 +      atype3 = atomTypeCont_.find(at3key);
180  
181 <  
181 >      // query atom types for their chains of responsibility
182 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
183 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
184 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
185  
186 +      std::vector<AtomType*>::iterator i;
187 +      std::vector<AtomType*>::iterator j;
188 +      std::vector<AtomType*>::iterator k;
189  
190 < BendType* ForceField::getMatchingBendType(const string &at1, const string &at2,
191 <                                          const string &at3);
192 < TorsionType* ForceField::getMatchingTorsionType(const string &at1, const string &at2,
193 <                                                const string &at3, const string &at4);
190 >      int ii = 0;
191 >      int jj = 0;
192 >      int kk = 0;
193 >      int IKscore;
194  
195 < double ForceField::getRcutForAtomType(AtomType* at);
195 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
196  
197 +      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
198 +        ii = 0;
199 +        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
200 +          kk = 0;
201 +          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
202 +          
203 +            IKscore = ii + kk;
204  
205 < vector<vector<string> > generateWildcardSequence(const vector<string> atomTypes) {
206 <  
207 <   vector<vector<string> > results;
205 >            std::vector<std::string> myKeys;
206 >            myKeys.push_back((*i)->getName());
207 >            myKeys.push_back((*j)->getName());
208 >            myKeys.push_back((*k)->getName());
209  
210 <  
210 >            BendType* bendType = bendTypeCont_.find(myKeys);
211 >            if (bendType) {
212 >              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
213 >            }
214 >            kk++;
215 >          }
216 >          ii++;
217 >        }
218 >        jj++;
219 >      }
220 >      
221 >      if (foundBends.size() > 0) {
222 >        std::sort(foundBends.begin(), foundBends.end());
223 >        int jscore = foundBends[0].first;
224 >        int ikscore = foundBends[0].second;
225 >        std::vector<std::string> theKeys = foundBends[0].third;      
226 >        
227 >        BendType* bestType = bendTypeCont_.find(theKeys);  
228 >        return bestType;
229 >      } else {        
230 >        //if no exact match found, try wild card match
231 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
232 >      }
233 >    }
234 >  }
235  
236 +  TorsionType* ForceField::getTorsionType(const std::string &at1,
237 +                                          const std::string &at2,
238 +                                          const std::string &at3,
239 +                                          const std::string &at4) {
240 +    std::vector<std::string> keys;
241 +    keys.push_back(at1);
242 +    keys.push_back(at2);    
243 +    keys.push_back(at3);    
244 +    keys.push_back(at4);    
245  
80   vector<vector< string> > getAllWildcardPermutations(const vector<string> myAts) {
81    
82     int nStrings;
83     vector<string> oneResult;
84     vector<vector<string> > allResults;
246  
247 <     nStrings = myAts.size();
247 >    //try exact match first
248 >    TorsionType* torsionType = torsionTypeCont_.find(keys);
249 >    if (torsionType) {
250 >      return torsionType;
251 >    } else {
252  
253 <     if (nStrings == 1) {
254 <       oneResult.push_back(wildcardCharacter);
255 <       allResults.push_back(oneResult);
256 <       return allResults;
257 <     } else {
258 <      
259 <       for (i=0; i < nStrings; i++) {
260 <         oneResult = myAts;
261 <         replace(oneResult.begin(), oneResult.end(),
253 >      AtomType* atype1;
254 >      AtomType* atype2;
255 >      AtomType* atype3;
256 >      AtomType* atype4;
257 >      std::vector<std::string> at1key;
258 >      at1key.push_back(at1);
259 >      atype1 = atomTypeCont_.find(at1key);
260 >  
261 >      std::vector<std::string> at2key;
262 >      at2key.push_back(at2);
263 >      atype2 = atomTypeCont_.find(at2key);
264 >
265 >      std::vector<std::string> at3key;
266 >      at3key.push_back(at3);
267 >      atype3 = atomTypeCont_.find(at3key);
268 >
269 >      std::vector<std::string> at4key;
270 >      at4key.push_back(at4);
271 >      atype4 = atomTypeCont_.find(at4key);
272 >
273 >      // query atom types for their chains of responsibility
274 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
275 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
276 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
277 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
278 >
279 >      std::vector<AtomType*>::iterator i;
280 >      std::vector<AtomType*>::iterator j;
281 >      std::vector<AtomType*>::iterator k;
282 >      std::vector<AtomType*>::iterator l;
283 >
284 >      int ii = 0;
285 >      int jj = 0;
286 >      int kk = 0;
287 >      int ll = 0;
288 >      int ILscore;
289 >      int JKscore;
290 >
291 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
292 >
293 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
294 >        kk = 0;
295 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
296 >          ii = 0;      
297 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
298 >            ll = 0;
299 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
300 >          
301 >              ILscore = ii + ll;
302 >              JKscore = jj + kk;
303 >
304 >              std::vector<std::string> myKeys;
305 >              myKeys.push_back((*i)->getName());
306 >              myKeys.push_back((*j)->getName());
307 >              myKeys.push_back((*k)->getName());
308 >              myKeys.push_back((*l)->getName());
309 >
310 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
311 >              if (torsionType) {
312 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
313 >              }
314 >              ll++;
315 >            }
316 >            ii++;
317 >          }
318 >          kk++;
319 >        }
320 >        jj++;
321 >      }
322 >      
323 >      if (foundTorsions.size() > 0) {
324 >        std::sort(foundTorsions.begin(), foundTorsions.end());
325 >        int jkscore = foundTorsions[0].first;
326 >        int ilscore = foundTorsions[0].second;
327 >        std::vector<std::string> theKeys = foundTorsions[0].third;
328 >        
329 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
330 >        return bestType;
331 >      } else {
332 >        //if no exact match found, try wild card match
333 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
334 >      }
335 >    }
336 >  }
337 >
338 >  InversionType* ForceField::getInversionType(const std::string &at1,
339 >                                              const std::string &at2,
340 >                                              const std::string &at3,
341 >                                              const std::string &at4) {
342 >    std::vector<std::string> keys;
343 >    keys.push_back(at1);
344 >    keys.push_back(at2);    
345 >    keys.push_back(at3);    
346 >    keys.push_back(at4);    
347 >
348 >    //try exact match first
349 >    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
350 >    if (inversionType) {
351 >      return inversionType;
352 >    } else {
353 >      
354 >      AtomType* atype1;
355 >      AtomType* atype2;
356 >      AtomType* atype3;
357 >      AtomType* atype4;
358 >      std::vector<std::string> at1key;
359 >      at1key.push_back(at1);
360 >      atype1 = atomTypeCont_.find(at1key);
361 >      
362 >      std::vector<std::string> at2key;
363 >      at2key.push_back(at2);
364 >      atype2 = atomTypeCont_.find(at2key);
365 >      
366 >      std::vector<std::string> at3key;
367 >      at3key.push_back(at3);
368 >      atype3 = atomTypeCont_.find(at3key);
369 >      
370 >      std::vector<std::string> at4key;
371 >      at4key.push_back(at4);
372 >      atype4 = atomTypeCont_.find(at4key);
373 >
374 >      // query atom types for their chains of responsibility
375 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
376 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
377 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
378 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
379 >
380 >      std::vector<AtomType*>::iterator i;
381 >      std::vector<AtomType*>::iterator j;
382 >      std::vector<AtomType*>::iterator k;
383 >      std::vector<AtomType*>::iterator l;
384 >
385 >      int ii = 0;
386 >      int jj = 0;
387 >      int kk = 0;
388 >      int ll = 0;
389 >      int Iscore;
390 >      int JKLscore;
391 >      
392 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
393 >      
394 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
395 >        kk = 0;
396 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
397 >          ii = 0;      
398 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
399 >            ll = 0;
400 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
401 >              
402 >              Iscore = ii;
403 >              JKLscore = jj + kk + ll;
404 >              
405 >              std::vector<std::string> myKeys;
406 >              myKeys.push_back((*i)->getName());
407 >              myKeys.push_back((*j)->getName());
408 >              myKeys.push_back((*k)->getName());
409 >              myKeys.push_back((*l)->getName());
410 >              
411 >              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
412 >              if (inversionType) {
413 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
414 >              }
415 >              ll++;
416 >            }
417 >            ii++;
418 >          }
419 >          kk++;
420 >        }
421 >        jj++;
422 >      }
423 >        
424 >      if (foundInversions.size() > 0) {
425 >        std::sort(foundInversions.begin(), foundInversions.end());
426 >        int iscore = foundInversions[0].first;
427 >        int jklscore = foundInversions[0].second;
428 >        std::vector<std::string> theKeys = foundInversions[0].third;
429 >        
430 >        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
431 >        return bestType;
432 >      } else {
433 >        //if no exact match found, try wild card match
434 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
435 >      }
436 >    }
437 >  }
438 >  
439 >  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
440 >    std::vector<std::string> keys;
441 >    keys.push_back(at1);
442 >    keys.push_back(at2);    
443 >    
444 >    //try exact match first
445 >    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
446 >    if (nbiType) {
447 >      return nbiType;
448 >    } else {
449 >      //if no exact match found, try wild card match
450 >      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
451 >    }    
452 >  }
453 >  
454 >  BondType* ForceField::getExactBondType(const std::string &at1,
455 >                                         const std::string &at2){
456 >    std::vector<std::string> keys;
457 >    keys.push_back(at1);
458 >    keys.push_back(at2);    
459 >    return bondTypeCont_.find(keys);
460 >  }
461 >  
462 >  BendType* ForceField::getExactBendType(const std::string &at1,
463 >                                         const std::string &at2,
464 >                                         const std::string &at3){
465 >    std::vector<std::string> keys;
466 >    keys.push_back(at1);
467 >    keys.push_back(at2);    
468 >    keys.push_back(at3);    
469 >    return bendTypeCont_.find(keys);
470 >  }
471 >  
472 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
473 >                                               const std::string &at2,
474 >                                               const std::string &at3,
475 >                                               const std::string &at4){
476 >    std::vector<std::string> keys;
477 >    keys.push_back(at1);
478 >    keys.push_back(at2);    
479 >    keys.push_back(at3);    
480 >    keys.push_back(at4);  
481 >    return torsionTypeCont_.find(keys);
482 >  }
483 >  
484 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
485 >                                                   const std::string &at2,
486 >                                                   const std::string &at3,
487 >                                                   const std::string &at4){
488 >    std::vector<std::string> keys;
489 >    keys.push_back(at1);
490 >    keys.push_back(at2);    
491 >    keys.push_back(at3);    
492 >    keys.push_back(at4);  
493 >    return inversionTypeCont_.find(keys);
494 >  }
495 >  
496 >  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
497 >    std::vector<std::string> keys;
498 >    keys.push_back(at1);
499 >    keys.push_back(at2);    
500 >    return nonBondedInteractionTypeCont_.find(keys);
501 >  }
502 >  
503 >
504 >  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
505 >    std::vector<std::string> keys;
506 >    keys.push_back(at);
507 >    return atomTypeCont_.add(keys, atomType);
508 >  }
509 >
510 >  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
511 >    std::vector<std::string> keys;
512 >    keys.push_back(at);
513 >    return atomTypeCont_.replace(keys, atomType);
514 >  }
515 >
516 >  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
517 >                               BondType* bondType) {
518 >    std::vector<std::string> keys;
519 >    keys.push_back(at1);
520 >    keys.push_back(at2);    
521 >    return bondTypeCont_.add(keys, bondType);    
522 >  }
523 >  
524 >  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
525 >                               const std::string &at3, BendType* bendType) {
526 >    std::vector<std::string> keys;
527 >    keys.push_back(at1);
528 >    keys.push_back(at2);    
529 >    keys.push_back(at3);    
530 >    return bendTypeCont_.add(keys, bendType);
531 >  }
532 >  
533 >  bool ForceField::addTorsionType(const std::string &at1,
534 >                                  const std::string &at2,
535 >                                  const std::string &at3,
536 >                                  const std::string &at4,
537 >                                  TorsionType* torsionType) {
538 >    std::vector<std::string> keys;
539 >    keys.push_back(at1);
540 >    keys.push_back(at2);    
541 >    keys.push_back(at3);    
542 >    keys.push_back(at4);    
543 >    return torsionTypeCont_.add(keys, torsionType);
544 >  }
545 >
546 >  bool ForceField::addInversionType(const std::string &at1,
547 >                                    const std::string &at2,
548 >                                    const std::string &at3,
549 >                                    const std::string &at4,
550 >                                    InversionType* inversionType) {
551 >    std::vector<std::string> keys;
552 >    keys.push_back(at1);
553 >    keys.push_back(at2);    
554 >    keys.push_back(at3);    
555 >    keys.push_back(at4);    
556 >    return inversionTypeCont_.add(keys, inversionType);
557 >  }
558 >  
559 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
560 >                                               const std::string &at2,
561 >                                               NonBondedInteractionType* nbiType) {
562 >    std::vector<std::string> keys;
563 >    keys.push_back(at1);
564 >    keys.push_back(at2);    
565 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
566 >  }
567 >  
568 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
569 >    /**@todo */
570 >    GenericData* data;
571 >    RealType rcut = 0.0;
572 >    
573 >    if (at->isLennardJones()) {
574 >      data = at->getPropertyByName("LennardJones");
575 >      if (data != NULL) {
576 >        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
577 >        
578 >        if (ljData != NULL) {
579 >          LJParam ljParam = ljData->getData();
580 >          
581 >          //by default use 2.5*sigma as cutoff radius
582 >          rcut = 2.5 * ljParam.sigma;
583 >          
584 >        } else {
585 >          sprintf( painCave.errMsg,
586 >                   "Can not cast GenericData to LJParam\n");
587 >          painCave.severity = OPENMD_ERROR;
588 >          painCave.isFatal = 1;
589 >          simError();          
590 >        }            
591 >      } else {
592 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
593 >        painCave.severity = OPENMD_ERROR;
594 >        painCave.isFatal = 1;
595 >        simError();          
596 >      }
597 >    }
598 >    return rcut;    
599 >  }
600 >  
601 >
602 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
603 >    std::string forceFieldFilename(filename);
604 >    ifstrstream* ffStream = new ifstrstream();
605 >    
606 >    //try to open the force filed file in current directory first    
607 >    ffStream->open(forceFieldFilename.c_str());
608 >    if(!ffStream->is_open()){
609 >
610 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
611 >      ffStream->open( forceFieldFilename.c_str() );
612 >
613 >      //if current directory does not contain the force field file,
614 >      //try to open it in the path        
615 >      if(!ffStream->is_open()){
616 >
617 >        sprintf( painCave.errMsg,
618 >                 "Error opening the force field parameter file:\n"
619 >                 "\t%s\n"
620 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
621 >                 "variable?\n",
622 >                 forceFieldFilename.c_str() );
623 >        painCave.severity = OPENMD_ERROR;
624 >        painCave.isFatal = 1;
625 >        simError();
626 >      }
627 >    }  
628 >    return ffStream;
629 >  }
630 >
631 > } //end namespace OpenMD

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 206 by gezelter, Thu Nov 4 20:51:23 2004 UTC vs.
branches/development/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 1530 by gezelter, Tue Dec 28 21:47:55 2010 UTC

# Line 0 | Line 1
1 + Author Id Revision Date

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines