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

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines