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

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines