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

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines