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

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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines