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

Comparing:
trunk/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 475 by tim, Tue Apr 12 18:30:37 2005 UTC vs.
branches/development/src/UseTheForce/ForceField.cpp (property svn:keywords), Revision 1532 by gezelter, Wed Dec 29 19:59:21 2010 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines