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 475 by tim, Tue Apr 12 18:30:37 2005 UTC vs.
Revision 1303 by cli2, Mon Oct 13 21:35:22 2008 UTC

# Line 1 | Line 1
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
# Line 39 | Line 39
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 <  */
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 < ForceField::ForceField() {
59 >  ForceField::ForceField() {
60      char* tempPath;
61      tempPath = getenv("FORCE_PARAM_PATH");
62  
63      if (tempPath == NULL) {
64 <        //convert a macro from compiler to a string in c++
65 <        STR_DEFINE(ffPath_, FRC_PATH );
64 >      //convert a macro from compiler to a string in c++
65 >      STR_DEFINE(ffPath_, FRC_PATH );
66      } else {
67 <        ffPath_ = tempPath;
67 >      ffPath_ = tempPath;
68      }
69 < }
69 >  }
70  
71  
72 < ForceField::~ForceField() {
72 >  ForceField::~ForceField() {
73      deleteAtypes();
74 < }
74 >    deleteSwitch();
75 >  }
76  
77 < AtomType* ForceField::getAtomType(const std::string &at) {
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 < }
81 >  }
82  
83 < BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
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);    
# Line 83 | Line 89 | BondType* ForceField::getBondType(const std::string &a
89      //try exact match first
90      BondType* bondType = bondTypeCont_.find(keys);
91      if (bondType) {
92 <        return bondType;
92 >      return bondType;
93      } else {
94 <        //if no exact match found, try wild card match
95 <        return bondTypeCont_.find(keys, wildCardAtomTypeName_);
96 <    }
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 < }
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 < BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
109 <                        const std::string &at3) {
108 >      std::vector<AtomType*>::iterator i;
109 >      std::vector<AtomType*>::iterator j;
110 >
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 >
137 >      if (foundBonds.size() > 0) {
138 >        // sort the foundBonds by the score:
139 >        std::sort(foundBonds.begin(), foundBonds.end());
140 >    
141 >        int bestScore = foundBonds[0].first;
142 >        std::vector<std::string> theKeys = foundBonds[0].second;
143 >        
144 >        BondType* bestType = bondTypeCont_.find(theKeys);
145 >        
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 >  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);    
# Line 101 | Line 162 | BendType* ForceField::getBendType(const std::string &a
162      //try exact match first
163      BendType* bendType = bendTypeCont_.find(keys);
164      if (bendType) {
165 <        return bendType;
165 >      return bendType;
166      } else {
106 –        //if no exact match found, try wild card match
107 –        return bendTypeCont_.find(keys, wildCardAtomTypeName_);
108 –    }
109 – }
167  
168 < TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
169 <                              const std::string &at3, const std::string &at4) {
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 >      // 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 >      int ii = 0;
193 >      int jj = 0;
194 >      int kk = 0;
195 >      int IKscore;
196 >
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 >            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 >            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 >      if (foundBends.size() > 0) {
224 >        std::sort(foundBends.begin(), foundBends.end());
225 >        int jscore = foundBends[0].first;
226 >        int ikscore = foundBends[0].second;
227 >        std::vector<std::string> theKeys = foundBends[0].third;      
228 >        
229 >        BendType* bestType = bendTypeCont_.find(theKeys);  
230 >        return bestType;
231 >      } else {        
232 >        //if no exact match found, try wild card match
233 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
234 >      }
235 >    }
236 >  }
237 >
238 >  TorsionType* ForceField::getTorsionType(const std::string &at1,
239 >                                          const std::string &at2,
240 >                                          const std::string &at3,
241 >                                          const std::string &at4) {
242      std::vector<std::string> keys;
243      keys.push_back(at1);
244      keys.push_back(at2);    
245      keys.push_back(at3);    
246      keys.push_back(at4);    
247  
248 +
249 +    //try exact match first
250      TorsionType* torsionType = torsionTypeCont_.find(keys);
251      if (torsionType) {
252 <        return torsionType;
252 >      return torsionType;
253      } else {
123 –        //if no exact match found, try wild card match
124 –        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
125 –    }
126 –    
127 –    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
254  
255 < }
255 >      AtomType* atype1;
256 >      AtomType* atype2;
257 >      AtomType* atype3;
258 >      AtomType* atype4;
259 >      std::vector<std::string> at1key;
260 >      at1key.push_back(at1);
261 >      atype1 = atomTypeCont_.find(at1key);
262 >  
263 >      std::vector<std::string> at2key;
264 >      at2key.push_back(at2);
265 >      atype2 = atomTypeCont_.find(at2key);
266  
267 < BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
267 >      std::vector<std::string> at3key;
268 >      at3key.push_back(at3);
269 >      atype3 = atomTypeCont_.find(at3key);
270 >
271 >      std::vector<std::string> at4key;
272 >      at4key.push_back(at4);
273 >      atype4 = atomTypeCont_.find(at4key);
274 >
275 >      // query atom types for their chains of responsibility
276 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
277 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
278 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
279 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
280 >
281 >      std::vector<AtomType*>::iterator i;
282 >      std::vector<AtomType*>::iterator j;
283 >      std::vector<AtomType*>::iterator k;
284 >      std::vector<AtomType*>::iterator l;
285 >
286 >      int ii = 0;
287 >      int jj = 0;
288 >      int kk = 0;
289 >      int ll = 0;
290 >      int ILscore;
291 >      int JKscore;
292 >
293 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
294 >
295 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
296 >        kk = 0;
297 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
298 >          ii = 0;      
299 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
300 >            ll = 0;
301 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
302 >          
303 >              ILscore = ii + ll;
304 >              JKscore = jj + kk;
305 >
306 >              std::vector<std::string> myKeys;
307 >              myKeys.push_back((*i)->getName());
308 >              myKeys.push_back((*j)->getName());
309 >              myKeys.push_back((*k)->getName());
310 >              myKeys.push_back((*l)->getName());
311 >
312 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
313 >              if (torsionType) {
314 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
315 >              }
316 >              ll++;
317 >            }
318 >            ii++;
319 >          }
320 >          kk++;
321 >        }
322 >        jj++;
323 >      }
324 >      
325 >      if (foundTorsions.size() > 0) {
326 >        std::sort(foundTorsions.begin(), foundTorsions.end());
327 >        int jkscore = foundTorsions[0].first;
328 >        int ilscore = foundTorsions[0].second;
329 >        std::vector<std::string> theKeys = foundTorsions[0].third;
330 >        
331 >        TorsionType* bestType = torsionTypeCont_.find(theKeys);
332 >        return bestType;
333 >      } else {
334 >        //if no exact match found, try wild card match
335 >        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
336 >      }
337 >    }
338 >  }
339 >
340 >  InversionType* ForceField::getInversionType(const std::string &at1,
341 >                                              const std::string &at2,
342 >                                              const std::string &at3,
343 >                                              const std::string &at4) {
344      std::vector<std::string> keys;
345      keys.push_back(at1);
346      keys.push_back(at2);    
347 <    return bondTypeCont_.find(keys);
348 < }
347 >    keys.push_back(at3);    
348 >    keys.push_back(at4);    
349  
350 < BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
351 <                            const std::string &at3){
350 >    //try exact match first
351 >    InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(keys);
352 >    if (inversionType) {
353 >      return inversionType;
354 >    } else {
355 >      
356 >      AtomType* atype1;
357 >      AtomType* atype2;
358 >      AtomType* atype3;
359 >      AtomType* atype4;
360 >      std::vector<std::string> at1key;
361 >      at1key.push_back(at1);
362 >      atype1 = atomTypeCont_.find(at1key);
363 >      
364 >      std::vector<std::string> at2key;
365 >      at2key.push_back(at2);
366 >      atype2 = atomTypeCont_.find(at2key);
367 >      
368 >      std::vector<std::string> at3key;
369 >      at3key.push_back(at3);
370 >      atype3 = atomTypeCont_.find(at3key);
371 >      
372 >      std::vector<std::string> at4key;
373 >      at4key.push_back(at4);
374 >      atype4 = atomTypeCont_.find(at4key);
375 >
376 >      // query atom types for their chains of responsibility
377 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
378 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
379 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
380 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
381 >
382 >      std::vector<AtomType*>::iterator i;
383 >      std::vector<AtomType*>::iterator j;
384 >      std::vector<AtomType*>::iterator k;
385 >      std::vector<AtomType*>::iterator l;
386 >
387 >      int ii = 0;
388 >      int jj = 0;
389 >      int kk = 0;
390 >      int ll = 0;
391 >      int Iscore;
392 >      int JKLscore;
393 >      
394 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
395 >      
396 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
397 >        kk = 0;
398 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
399 >          ii = 0;      
400 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
401 >            ll = 0;
402 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
403 >              
404 >              Iscore = ii;
405 >              JKLscore = jj + kk + ll;
406 >              
407 >              std::vector<std::string> myKeys;
408 >              myKeys.push_back((*i)->getName());
409 >              myKeys.push_back((*j)->getName());
410 >              myKeys.push_back((*k)->getName());
411 >              myKeys.push_back((*l)->getName());
412 >              
413 >              InversionType* inversionType = inversionTypeCont_.permutedFindSkippingFirstElement(myKeys);
414 >              if (inversionType) {
415 >                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
416 >              }
417 >              ll++;
418 >            }
419 >            ii++;
420 >          }
421 >          kk++;
422 >        }
423 >        jj++;
424 >      }
425 >        
426 >      if (foundInversions.size() > 0) {
427 >        std::sort(foundInversions.begin(), foundInversions.end());
428 >        int iscore = foundInversions[0].first;
429 >        int jklscore = foundInversions[0].second;
430 >        std::vector<std::string> theKeys = foundInversions[0].third;
431 >        
432 >        InversionType* bestType = inversionTypeCont_.permutedFindSkippingFirstElement(theKeys);
433 >        return bestType;
434 >      } else {
435 >        //if no exact match found, try wild card match
436 >        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
437 >      }
438 >    }
439 >  }
440 >  
441 >  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
442      std::vector<std::string> keys;
443      keys.push_back(at1);
444      keys.push_back(at2);    
445 +    
446 +    //try exact match first
447 +    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
448 +    if (nbiType) {
449 +      return nbiType;
450 +    } else {
451 +      //if no exact match found, try wild card match
452 +      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
453 +    }    
454 +  }
455 +  
456 +  BondType* ForceField::getExactBondType(const std::string &at1,
457 +                                         const std::string &at2){
458 +    std::vector<std::string> keys;
459 +    keys.push_back(at1);
460 +    keys.push_back(at2);    
461 +    return bondTypeCont_.find(keys);
462 +  }
463 +  
464 +  BendType* ForceField::getExactBendType(const std::string &at1,
465 +                                         const std::string &at2,
466 +                                         const std::string &at3){
467 +    std::vector<std::string> keys;
468 +    keys.push_back(at1);
469 +    keys.push_back(at2);    
470      keys.push_back(at3);    
471      return bendTypeCont_.find(keys);
472 < }
473 <
474 < TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
475 <                                  const std::string &at3, const std::string &at4){
472 >  }
473 >  
474 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
475 >                                               const std::string &at2,
476 >                                               const std::string &at3,
477 >                                               const std::string &at4){
478      std::vector<std::string> keys;
479      keys.push_back(at1);
480      keys.push_back(at2);    
481      keys.push_back(at3);    
482      keys.push_back(at4);  
483      return torsionTypeCont_.find(keys);
484 < }
485 < bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
484 >  }
485 >  
486 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
487 >                                                   const std::string &at2,
488 >                                                   const std::string &at3,
489 >                                                   const std::string &at4){
490      std::vector<std::string> keys;
491 +    keys.push_back(at1);
492 +    keys.push_back(at2);    
493 +    keys.push_back(at3);    
494 +    keys.push_back(at4);  
495 +    return inversionTypeCont_.find(keys);
496 +  }
497 +  
498 +  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
499 +    std::vector<std::string> keys;
500 +    keys.push_back(at1);
501 +    keys.push_back(at2);    
502 +    return nonBondedInteractionTypeCont_.find(keys);
503 +  }
504 +  
505 +
506 +  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
507 +    std::vector<std::string> keys;
508      keys.push_back(at);
509      return atomTypeCont_.add(keys, atomType);
510 < }
510 >  }
511  
512 < bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
512 >  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
513      std::vector<std::string> keys;
514 +    keys.push_back(at);
515 +    return atomTypeCont_.replace(keys, atomType);
516 +  }
517 +
518 +  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
519 +                               BondType* bondType) {
520 +    std::vector<std::string> keys;
521      keys.push_back(at1);
522      keys.push_back(at2);    
523 <    return bondTypeCont_.add(keys, bondType);
524 <
525 < }
526 <
527 < bool ForceField::addBendType(const std::string &at1, const std::string &at2,
171 <                        const std::string &at3, BendType* bendType) {
523 >    return bondTypeCont_.add(keys, bondType);    
524 >  }
525 >  
526 >  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
527 >                               const std::string &at3, BendType* bendType) {
528      std::vector<std::string> keys;
529      keys.push_back(at1);
530      keys.push_back(at2);    
531      keys.push_back(at3);    
532      return bendTypeCont_.add(keys, bendType);
533 < }
534 <
535 < bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
536 <                              const std::string &at3, const std::string &at4, TorsionType* torsionType) {
533 >  }
534 >  
535 >  bool ForceField::addTorsionType(const std::string &at1,
536 >                                  const std::string &at2,
537 >                                  const std::string &at3,
538 >                                  const std::string &at4,
539 >                                  TorsionType* torsionType) {
540      std::vector<std::string> keys;
541      keys.push_back(at1);
542      keys.push_back(at2);    
543      keys.push_back(at3);    
544      keys.push_back(at4);    
545      return torsionTypeCont_.add(keys, torsionType);
546 < }
546 >  }
547  
548 < double ForceField::getRcutFromAtomType(AtomType* at) {
548 >  bool ForceField::addInversionType(const std::string &at1,
549 >                                    const std::string &at2,
550 >                                    const std::string &at3,
551 >                                    const std::string &at4,
552 >                                    InversionType* inversionType) {
553 >    std::vector<std::string> keys;
554 >    keys.push_back(at1);
555 >    keys.push_back(at2);    
556 >    keys.push_back(at3);    
557 >    keys.push_back(at4);    
558 >    return inversionTypeCont_.add(keys, inversionType);
559 >  }
560 >  
561 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
562 >                                               const std::string &at2,
563 >                                               NonBondedInteractionType* nbiType) {
564 >    std::vector<std::string> keys;
565 >    keys.push_back(at1);
566 >    keys.push_back(at2);    
567 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
568 >  }
569 >  
570 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
571      /**@todo */
572      GenericData* data;
573 <    double rcut = 0.0;
574 <
573 >    RealType rcut = 0.0;
574 >    
575      if (at->isLennardJones()) {
576 <        data = at->getPropertyByName("LennardJones");
577 <        if (data != NULL) {
578 <            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
579 <
580 <            if (ljData != NULL) {
581 <                LJParam ljParam = ljData->getData();
582 <
583 <                //by default use 2.5*sigma as cutoff radius
584 <                rcut = 2.5 * ljParam.sigma;
585 <                
586 <            } else {
587 <                    sprintf( painCave.errMsg,
588 <                           "Can not cast GenericData to LJParam\n");
589 <                    painCave.severity = OOPSE_ERROR;
590 <                    painCave.isFatal = 1;
591 <                    simError();          
592 <            }            
593 <        } else {
594 <            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
595 <            painCave.severity = OOPSE_ERROR;
596 <            painCave.isFatal = 1;
597 <            simError();          
598 <        }
576 >      data = at->getPropertyByName("LennardJones");
577 >      if (data != NULL) {
578 >        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
579 >        
580 >        if (ljData != NULL) {
581 >          LJParam ljParam = ljData->getData();
582 >          
583 >          //by default use 2.5*sigma as cutoff radius
584 >          rcut = 2.5 * ljParam.sigma;
585 >          
586 >        } else {
587 >          sprintf( painCave.errMsg,
588 >                   "Can not cast GenericData to LJParam\n");
589 >          painCave.severity = OOPSE_ERROR;
590 >          painCave.isFatal = 1;
591 >          simError();          
592 >        }            
593 >      } else {
594 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
595 >        painCave.severity = OOPSE_ERROR;
596 >        painCave.isFatal = 1;
597 >        simError();          
598 >      }
599      }
219 –
600      return rcut;    
601 < }
601 >  }
602 >  
603  
604 <
224 < ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
604 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
605      std::string forceFieldFilename(filename);
606      ifstrstream* ffStream = new ifstrstream();
607      
# Line 229 | Line 609 | ifstrstream* ForceField::openForceFieldFile(const std:
609      ffStream->open(forceFieldFilename.c_str());
610      if(!ffStream->is_open()){
611  
612 <        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
613 <        ffStream->open( forceFieldFilename.c_str() );
612 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
613 >      ffStream->open( forceFieldFilename.c_str() );
614  
615 <        //if current directory does not contain the force field file,
616 <        //try to open it in the path        
617 <        if(!ffStream->is_open()){
615 >      //if current directory does not contain the force field file,
616 >      //try to open it in the path        
617 >      if(!ffStream->is_open()){
618  
619 <            sprintf( painCave.errMsg,
620 <               "Error opening the force field parameter file:\n"
621 <               "\t%s\n"
622 <               "\tHave you tried setting the FORCE_PARAM_PATH environment "
623 <               "variable?\n",
624 <               forceFieldFilename.c_str() );
625 <            painCave.severity = OOPSE_ERROR;
626 <            painCave.isFatal = 1;
627 <            simError();
628 <        }
619 >        sprintf( painCave.errMsg,
620 >                 "Error opening the force field parameter file:\n"
621 >                 "\t%s\n"
622 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
623 >                 "variable?\n",
624 >                 forceFieldFilename.c_str() );
625 >        painCave.severity = OOPSE_ERROR;
626 >        painCave.isFatal = 1;
627 >        simError();
628 >      }
629      }  
250 –
630      return ffStream;
631 +  }
632  
633 < }
634 <
633 >  void ForceField::setFortranForceOptions(){
634 >    ForceOptions theseFortranOptions;
635 >    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
636 >    setfForceOptions(&theseFortranOptions);
637 >  }
638   } //end namespace oopse

Diff Legend

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