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

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines