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 1629 by gezelter, Wed Sep 14 21:15:17 2011 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 {
106        //if no exact match found, try wild card match
107        return bendTypeCont_.find(keys, wildCardAtomTypeName_);
108    }
109 }
176  
177 < TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
178 <                              const std::string &at3, const std::string &at4) {
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 >  }
246 >
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 >    
452      std::vector<std::string> keys;
453      keys.push_back(at1);
454      keys.push_back(at2);    
455 +    
456 +    //try exact match first
457 +    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
458 +    if (nbiType) {
459 +      return nbiType;
460 +    } else {
461 +      AtomType* atype1;
462 +      AtomType* atype2;
463 +      std::vector<std::string> at1key;
464 +      at1key.push_back(at1);
465 +      atype1 = atomTypeCont_.find(at1key);
466 +      
467 +      std::vector<std::string> at2key;
468 +      at2key.push_back(at2);
469 +      atype2 = atomTypeCont_.find(at2key);
470 +      
471 +      // query atom types for their chains of responsibility
472 +      std::vector<AtomType*> at1Chain = atype1->allYourBase();
473 +      std::vector<AtomType*> at2Chain = atype2->allYourBase();
474 +      
475 +      std::vector<AtomType*>::iterator i;
476 +      std::vector<AtomType*>::iterator j;
477 +      
478 +      int ii = 0;
479 +      int jj = 0;
480 +      int nbiTypeScore;
481 +      
482 +      std::vector<std::pair<int, std::vector<std::string> > > foundNBI;
483 +      
484 +      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
485 +        jj = 0;
486 +        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
487 +          
488 +          nbiTypeScore = ii + jj;
489 +          
490 +          std::vector<std::string> myKeys;
491 +          myKeys.push_back((*i)->getName());
492 +          myKeys.push_back((*j)->getName());
493 +          
494 +          NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(myKeys);
495 +          if (nbiType) {
496 +            foundNBI.push_back(std::make_pair(nbiTypeScore, myKeys));
497 +          }
498 +          jj++;
499 +        }
500 +        ii++;
501 +      }
502 +      
503 +      
504 +      if (foundNBI.size() > 0) {
505 +        // sort the foundNBI by the score:
506 +        std::sort(foundNBI.begin(), foundNBI.end());
507 +        
508 +        int bestScore = foundNBI[0].first;
509 +        std::vector<std::string> theKeys = foundNBI[0].second;
510 +        
511 +        NonBondedInteractionType* bestType = nonBondedInteractionTypeCont_.find(theKeys);        
512 +        return bestType;
513 +      } else {
514 +        //if no exact match found, try wild card match
515 +        return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
516 +      }
517 +    }
518 +  }
519 +  
520 +  BondType* ForceField::getExactBondType(const std::string &at1,
521 +                                         const std::string &at2){
522 +    std::vector<std::string> keys;
523 +    keys.push_back(at1);
524 +    keys.push_back(at2);    
525 +    return bondTypeCont_.find(keys);
526 +  }
527 +  
528 +  BendType* ForceField::getExactBendType(const std::string &at1,
529 +                                         const std::string &at2,
530 +                                         const std::string &at3){
531 +    std::vector<std::string> keys;
532 +    keys.push_back(at1);
533 +    keys.push_back(at2);    
534      keys.push_back(at3);    
535      return bendTypeCont_.find(keys);
536 < }
537 <
538 < TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
539 <                                  const std::string &at3, const std::string &at4){
536 >  }
537 >  
538 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
539 >                                               const std::string &at2,
540 >                                               const std::string &at3,
541 >                                               const std::string &at4){
542      std::vector<std::string> keys;
543      keys.push_back(at1);
544      keys.push_back(at2);    
545      keys.push_back(at3);    
546      keys.push_back(at4);  
547      return torsionTypeCont_.find(keys);
548 < }
549 < bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
548 >  }
549 >  
550 >  InversionType* ForceField::getExactInversionType(const std::string &at1,
551 >                                                   const std::string &at2,
552 >                                                   const std::string &at3,
553 >                                                   const std::string &at4){
554      std::vector<std::string> keys;
555 +    keys.push_back(at1);
556 +    keys.push_back(at2);    
557 +    keys.push_back(at3);    
558 +    keys.push_back(at4);  
559 +    return inversionTypeCont_.find(keys);
560 +  }
561 +  
562 +  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
563 +    std::vector<std::string> keys;
564 +    keys.push_back(at1);
565 +    keys.push_back(at2);    
566 +    return nonBondedInteractionTypeCont_.find(keys);
567 +  }
568 +  
569 +
570 +  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
571 +    std::vector<std::string> keys;
572      keys.push_back(at);
573 +    atypeIdentToName[atomType->getIdent()] = at;
574      return atomTypeCont_.add(keys, atomType);
575 < }
575 >  }
576  
577 < bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
577 >  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
578      std::vector<std::string> keys;
579 +    keys.push_back(at);
580 +    atypeIdentToName[atomType->getIdent()] = at;
581 +    return atomTypeCont_.replace(keys, atomType);
582 +  }
583 +
584 +  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
585 +                               BondType* bondType) {
586 +    std::vector<std::string> keys;
587      keys.push_back(at1);
588      keys.push_back(at2);    
589 <    return bondTypeCont_.add(keys, bondType);
590 <
591 < }
592 <
593 < bool ForceField::addBendType(const std::string &at1, const std::string &at2,
171 <                        const std::string &at3, BendType* bendType) {
589 >    return bondTypeCont_.add(keys, bondType);    
590 >  }
591 >  
592 >  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
593 >                               const std::string &at3, BendType* bendType) {
594      std::vector<std::string> keys;
595      keys.push_back(at1);
596      keys.push_back(at2);    
597      keys.push_back(at3);    
598      return bendTypeCont_.add(keys, bendType);
599 < }
600 <
601 < bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
602 <                              const std::string &at3, const std::string &at4, TorsionType* torsionType) {
599 >  }
600 >  
601 >  bool ForceField::addTorsionType(const std::string &at1,
602 >                                  const std::string &at2,
603 >                                  const std::string &at3,
604 >                                  const std::string &at4,
605 >                                  TorsionType* torsionType) {
606      std::vector<std::string> keys;
607      keys.push_back(at1);
608      keys.push_back(at2);    
609      keys.push_back(at3);    
610      keys.push_back(at4);    
611      return torsionTypeCont_.add(keys, torsionType);
612 < }
612 >  }
613  
614 < double ForceField::getRcutFromAtomType(AtomType* at) {
614 >  bool ForceField::addInversionType(const std::string &at1,
615 >                                    const std::string &at2,
616 >                                    const std::string &at3,
617 >                                    const std::string &at4,
618 >                                    InversionType* inversionType) {
619 >    std::vector<std::string> keys;
620 >    keys.push_back(at1);
621 >    keys.push_back(at2);    
622 >    keys.push_back(at3);    
623 >    keys.push_back(at4);    
624 >    return inversionTypeCont_.add(keys, inversionType);
625 >  }
626 >  
627 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
628 >                                               const std::string &at2,
629 >                                               NonBondedInteractionType* nbiType) {
630 >    std::vector<std::string> keys;
631 >    keys.push_back(at1);
632 >    keys.push_back(at2);    
633 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
634 >  }
635 >  
636 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
637      /**@todo */
638      GenericData* data;
639 <    double rcut = 0.0;
640 <
639 >    RealType rcut = 0.0;
640 >    
641      if (at->isLennardJones()) {
642 <        data = at->getPropertyByName("LennardJones");
643 <        if (data != NULL) {
644 <            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
645 <
646 <            if (ljData != NULL) {
647 <                LJParam ljParam = ljData->getData();
648 <
649 <                //by default use 2.5*sigma as cutoff radius
650 <                rcut = 2.5 * ljParam.sigma;
651 <                
652 <            } else {
653 <                    sprintf( painCave.errMsg,
654 <                           "Can not cast GenericData to LJParam\n");
655 <                    painCave.severity = OOPSE_ERROR;
656 <                    painCave.isFatal = 1;
657 <                    simError();          
658 <            }            
659 <        } else {
660 <            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
661 <            painCave.severity = OOPSE_ERROR;
662 <            painCave.isFatal = 1;
663 <            simError();          
664 <        }
642 >      data = at->getPropertyByName("LennardJones");
643 >      if (data != NULL) {
644 >        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
645 >        
646 >        if (ljData != NULL) {
647 >          LJParam ljParam = ljData->getData();
648 >          
649 >          //by default use 2.5*sigma as cutoff radius
650 >          rcut = 2.5 * ljParam.sigma;
651 >          
652 >        } else {
653 >          sprintf( painCave.errMsg,
654 >                   "Can not cast GenericData to LJParam\n");
655 >          painCave.severity = OPENMD_ERROR;
656 >          painCave.isFatal = 1;
657 >          simError();          
658 >        }            
659 >      } else {
660 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
661 >        painCave.severity = OPENMD_ERROR;
662 >        painCave.isFatal = 1;
663 >        simError();          
664 >      }
665      }
219
666      return rcut;    
667 < }
667 >  }
668 >  
669  
670 <
224 < ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
670 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
671      std::string forceFieldFilename(filename);
672      ifstrstream* ffStream = new ifstrstream();
673      
# Line 229 | Line 675 | ifstrstream* ForceField::openForceFieldFile(const std:
675      ffStream->open(forceFieldFilename.c_str());
676      if(!ffStream->is_open()){
677  
678 <        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
679 <        ffStream->open( forceFieldFilename.c_str() );
678 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
679 >      ffStream->open( forceFieldFilename.c_str() );
680  
681 <        //if current directory does not contain the force field file,
682 <        //try to open it in the path        
683 <        if(!ffStream->is_open()){
681 >      //if current directory does not contain the force field file,
682 >      //try to open it in the path        
683 >      if(!ffStream->is_open()){
684  
685 <            sprintf( painCave.errMsg,
686 <               "Error opening the force field parameter file:\n"
687 <               "\t%s\n"
688 <               "\tHave you tried setting the FORCE_PARAM_PATH environment "
689 <               "variable?\n",
690 <               forceFieldFilename.c_str() );
691 <            painCave.severity = OOPSE_ERROR;
692 <            painCave.isFatal = 1;
693 <            simError();
694 <        }
685 >        sprintf( painCave.errMsg,
686 >                 "Error opening the force field parameter file:\n"
687 >                 "\t%s\n"
688 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
689 >                 "variable?\n",
690 >                 forceFieldFilename.c_str() );
691 >        painCave.severity = OPENMD_ERROR;
692 >        painCave.isFatal = 1;
693 >        simError();
694 >      }
695      }  
250
696      return ffStream;
697 +  }
698  
699 < }
254 <
255 < } //end namespace oopse
699 > } //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 1629 by gezelter, Wed Sep 14 21:15:17 2011 UTC

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines