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 246 by gezelter, Wed Jan 12 22:41:40 2005 UTC vs.
Revision 1269 by gezelter, Tue Jul 1 13:28:23 2008 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 39 | Line 39
39   * such damages.
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 + #include "UseTheForce/DarkSide/fForceOptions_interface.h"
56 + #include "UseTheForce/DarkSide/switcheroo_interface.h"
57   namespace oopse {
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++
65 <        STR_DEFINE(ffPath_, FRC_PATH );
64 >      //convert a macro from compiler to a string in c++
65 >      STR_DEFINE(ffPath_, FRC_PATH );
66      } else {
67 <        ffPath_ = tempPath;
67 >      ffPath_ = tempPath;
68      }
69 < }
69 >  }
70  
71 < AtomType* ForceField::getAtomType(const std::string &at) {
71 >
72 >  ForceField::~ForceField() {
73 >    deleteAtypes();
74 >    deleteSwitch();
75 >  }
76 >
77 >  AtomType* ForceField::getAtomType(const std::string &at) {
78      std::vector<std::string> keys;
79      keys.push_back(at);
80      return atomTypeCont_.find(keys);
81 < }
81 >  }
82  
83 < BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
83 >  BondType* ForceField::getBondType(const std::string &at1,
84 >                                    const std::string &at2) {
85      std::vector<std::string> keys;
86      keys.push_back(at1);
87      keys.push_back(at2);    
# Line 77 | Line 89 | BondType* ForceField::getBondType(const std::string &a
89      //try exact match first
90      BondType* bondType = bondTypeCont_.find(keys);
91      if (bondType) {
92 <        return bondType;
92 >      return bondType;
93      } else {
94 <        //if no exact match found, try wild card match
95 <        return bondTypeCont_.find(keys, wildCardAtomTypeName_);
96 <    }
94 >      AtomType* atype1;
95 >      AtomType* atype2;
96 >      std::vector<std::string> at1key;
97 >      at1key.push_back(at1);
98 >      atype1 = atomTypeCont_.find(at1key);
99 >  
100 >      std::vector<std::string> at2key;
101 >      at2key.push_back(at2);
102 >      atype2 = atomTypeCont_.find(at2key);
103  
104 < }
104 >      // query atom types for their chains of responsibility
105 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
106 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
107  
108 < BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
109 <                        const std::string &at3) {
108 >      std::vector<AtomType*>::iterator i;
109 >      std::vector<AtomType*>::iterator j;
110 >
111 >      int ii = 0;
112 >      int jj = 0;
113 >      int bondTypeScore;
114 >
115 >      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
116 >
117 >      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
118 >        jj = 0;
119 >        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
120 >
121 >          bondTypeScore = ii + jj;
122 >
123 >          std::vector<std::string> myKeys;
124 >          myKeys.push_back((*i)->getName());
125 >          myKeys.push_back((*j)->getName());
126 >
127 >          BondType* bondType = bondTypeCont_.find(myKeys);
128 >          if (bondType) {
129 >            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
130 >          }
131 >          jj++;
132 >        }
133 >        ii++;
134 >      }
135 >
136 >      // sort the foundBonds by the score:
137 >
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 >      std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";      
144 >      BondType* bestType = bondTypeCont_.find(theKeys);
145 >      if (bestType)
146 >        return bestType;
147 >      else {
148 >        //if no exact match found, try wild card match
149 >        return bondTypeCont_.find(keys, wildCardAtomTypeName_);      
150 >      }
151 >    }
152 >  }
153 >  
154 >  BendType* ForceField::getBendType(const std::string &at1,
155 >                                    const std::string &at2,
156 >                                    const std::string &at3) {
157      std::vector<std::string> keys;
158      keys.push_back(at1);
159      keys.push_back(at2);    
# Line 95 | Line 162 | BendType* ForceField::getBendType(const std::string &a
162      //try exact match first
163      BendType* bendType = bendTypeCont_.find(keys);
164      if (bendType) {
165 <        return bendType;
165 >      return bendType;
166      } else {
167 <        //if no exact match found, try wild card match
168 <        return bendTypeCont_.find(keys, wildCardAtomTypeName_);
167 >
168 >      AtomType* atype1;
169 >      AtomType* atype2;
170 >      AtomType* atype3;
171 >      std::vector<std::string> at1key;
172 >      at1key.push_back(at1);
173 >      atype1 = atomTypeCont_.find(at1key);
174 >  
175 >      std::vector<std::string> at2key;
176 >      at2key.push_back(at2);
177 >      atype2 = atomTypeCont_.find(at2key);
178 >
179 >      std::vector<std::string> at3key;
180 >      at3key.push_back(at3);
181 >      atype3 = atomTypeCont_.find(at3key);
182 >
183 >      // query atom types for their chains of responsibility
184 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
185 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
186 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
187 >
188 >      std::vector<AtomType*>::iterator i;
189 >      std::vector<AtomType*>::iterator j;
190 >      std::vector<AtomType*>::iterator k;
191 >
192 >      int ii = 0;
193 >      int jj = 0;
194 >      int kk = 0;
195 >      int IKscore;
196 >
197 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
198 >
199 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
200 >        ii = 0;
201 >        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
202 >          kk = 0;
203 >          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
204 >          
205 >            IKscore = ii + kk;
206 >
207 >            std::vector<std::string> myKeys;
208 >            myKeys.push_back((*i)->getName());
209 >            myKeys.push_back((*j)->getName());
210 >            myKeys.push_back((*k)->getName());
211 >
212 >            BendType* bendType = bendTypeCont_.find(myKeys);
213 >            if (bendType) {
214 >              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
215 >            }
216 >            kk++;
217 >          }
218 >          ii++;
219 >        }
220 >        jj++;
221 >      }
222 >      
223 >      std::sort(foundBends.begin(), foundBends.end());
224 >
225 >      int jscore = foundBends[0].first;
226 >      int ikscore = foundBends[0].second;
227 >      std::vector<std::string> theKeys = foundBends[0].third;
228 >
229 >      std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";      
230 >
231 >      BendType* bestType = bendTypeCont_.find(theKeys);  
232 >      if (bestType)
233 >        return bestType;
234 >      else {
235 >      
236 >        //if no exact match found, try wild card match
237 >        return bendTypeCont_.find(keys, wildCardAtomTypeName_);      
238 >      }
239      }
240 < }
240 >  }
241  
242 < TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
243 <                              const std::string &at3, const std::string &at4) {
242 >
243 >  TorsionType* ForceField::getTorsionType(const std::string &at1,
244 >                                          const std::string &at2,
245 >                                          const std::string &at3,
246 >                                          const std::string &at4) {
247      std::vector<std::string> keys;
248      keys.push_back(at1);
249      keys.push_back(at2);    
250      keys.push_back(at3);    
251      keys.push_back(at4);    
252  
253 +
254 +    //try exact match first
255      TorsionType* torsionType = torsionTypeCont_.find(keys);
256      if (torsionType) {
257 <        return torsionType;
257 >      return torsionType;
258      } else {
117 –        //if no exact match found, try wild card match
118 –        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
119 –    }
120 –    
121 –    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
259  
260 < }
260 >      AtomType* atype1;
261 >      AtomType* atype2;
262 >      AtomType* atype3;
263 >      AtomType* atype4;
264 >      std::vector<std::string> at1key;
265 >      at1key.push_back(at1);
266 >      atype1 = atomTypeCont_.find(at1key);
267 >  
268 >      std::vector<std::string> at2key;
269 >      at2key.push_back(at2);
270 >      atype2 = atomTypeCont_.find(at2key);
271  
272 < BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
272 >      std::vector<std::string> at3key;
273 >      at3key.push_back(at3);
274 >      atype3 = atomTypeCont_.find(at3key);
275 >
276 >      std::vector<std::string> at4key;
277 >      at4key.push_back(at4);
278 >      atype4 = atomTypeCont_.find(at4key);
279 >
280 >      // query atom types for their chains of responsibility
281 >      std::vector<AtomType*> at1Chain = atype1->allYourBase();
282 >      std::vector<AtomType*> at2Chain = atype2->allYourBase();
283 >      std::vector<AtomType*> at3Chain = atype3->allYourBase();
284 >      std::vector<AtomType*> at4Chain = atype4->allYourBase();
285 >
286 >      std::vector<AtomType*>::iterator i;
287 >      std::vector<AtomType*>::iterator j;
288 >      std::vector<AtomType*>::iterator k;
289 >      std::vector<AtomType*>::iterator l;
290 >
291 >      int ii = 0;
292 >      int jj = 0;
293 >      int kk = 0;
294 >      int ll = 0;
295 >      int ILscore;
296 >      int JKscore;
297 >
298 >      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
299 >
300 >      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
301 >        kk = 0;
302 >        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
303 >          ii = 0;      
304 >          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
305 >            ll = 0;
306 >            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
307 >          
308 >              ILscore = ii + ll;
309 >              JKscore = jj + kk;
310 >
311 >              std::vector<std::string> myKeys;
312 >              myKeys.push_back((*i)->getName());
313 >              myKeys.push_back((*j)->getName());
314 >              myKeys.push_back((*k)->getName());
315 >              myKeys.push_back((*l)->getName());
316 >
317 >              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
318 >              if (torsionType) {
319 >                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
320 >              }
321 >              ll++;
322 >            }
323 >            ii++;
324 >          }
325 >          kk++;
326 >        }
327 >        jj++;
328 >      }
329 >      
330 >      std::sort(foundTorsions.begin(), foundTorsions.end());
331 >
332 >      int jkscore = foundTorsions[0].first;
333 >      int ilscore = foundTorsions[0].second;
334 >      std::vector<std::string> theKeys = foundTorsions[0].third;
335 >
336 >      std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
337 >
338 >      
339 >      TorsionType* bestType = torsionTypeCont_.find(theKeys);
340 >      if (bestType) {
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 >  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
350      std::vector<std::string> keys;
351      keys.push_back(at1);
352      keys.push_back(at2);    
353 +    
354 +    //try exact match first
355 +    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
356 +    if (nbiType) {
357 +      return nbiType;
358 +    } else {
359 +      //if no exact match found, try wild card match
360 +      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
361 +    }    
362 +  }
363 +  
364 +  BondType* ForceField::getExactBondType(const std::string &at1,
365 +                                         const std::string &at2){
366 +    std::vector<std::string> keys;
367 +    keys.push_back(at1);
368 +    keys.push_back(at2);    
369      return bondTypeCont_.find(keys);
370 < }
371 <
372 < BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
373 <                            const std::string &at3){
370 >  }
371 >  
372 >  BendType* ForceField::getExactBendType(const std::string &at1,
373 >                                         const std::string &at2,
374 >                                         const std::string &at3){
375      std::vector<std::string> keys;
376      keys.push_back(at1);
377      keys.push_back(at2);    
378      keys.push_back(at3);    
379      return bendTypeCont_.find(keys);
380 < }
381 <
382 < TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
383 <                                  const std::string &at3, const std::string &at4){
380 >  }
381 >  
382 >  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
383 >                                               const std::string &at2,
384 >                                               const std::string &at3,
385 >                                               const std::string &at4){
386      std::vector<std::string> keys;
387      keys.push_back(at1);
388      keys.push_back(at2);    
389      keys.push_back(at3);    
390      keys.push_back(at4);  
391      return torsionTypeCont_.find(keys);
392 < }
393 < bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
392 >  }
393 >
394 >  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
395      std::vector<std::string> keys;
396 +    keys.push_back(at1);
397 +    keys.push_back(at2);    
398 +    return nonBondedInteractionTypeCont_.find(keys);
399 +  }
400 +
401 +
402 +  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
403 +    std::vector<std::string> keys;
404      keys.push_back(at);
405      return atomTypeCont_.add(keys, atomType);
406 < }
406 >  }
407  
408 < bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
408 >  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
409 >                               BondType* bondType) {
410      std::vector<std::string> keys;
411      keys.push_back(at1);
412      keys.push_back(at2);    
413 <    return bondTypeCont_.add(keys, bondType);
414 <
415 < }
416 <
417 < bool ForceField::addBendType(const std::string &at1, const std::string &at2,
165 <                        const std::string &at3, BendType* bendType) {
413 >    return bondTypeCont_.add(keys, bondType);    
414 >  }
415 >  
416 >  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
417 >                               const std::string &at3, BendType* bendType) {
418      std::vector<std::string> keys;
419      keys.push_back(at1);
420      keys.push_back(at2);    
421      keys.push_back(at3);    
422      return bendTypeCont_.add(keys, bendType);
423 < }
424 <
425 < bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
426 <                              const std::string &at3, const std::string &at4, TorsionType* torsionType) {
423 >  }
424 >  
425 >  bool ForceField::addTorsionType(const std::string &at1,
426 >                                  const std::string &at2,
427 >                                  const std::string &at3,
428 >                                  const std::string &at4,
429 >                                  TorsionType* torsionType) {
430      std::vector<std::string> keys;
431      keys.push_back(at1);
432      keys.push_back(at2);    
433      keys.push_back(at3);    
434      keys.push_back(at4);    
435      return torsionTypeCont_.add(keys, torsionType);
436 < }
436 >  }
437  
438 < double ForceField::getRcutFromAtomType(AtomType* at) {
438 >  bool ForceField::addNonBondedInteractionType(const std::string &at1,
439 >                                               const std::string &at2,
440 >                                               NonBondedInteractionType* nbiType) {
441 >    std::vector<std::string> keys;
442 >    keys.push_back(at1);
443 >    keys.push_back(at2);    
444 >    return nonBondedInteractionTypeCont_.add(keys, nbiType);
445 >  }
446 >  
447 >  RealType ForceField::getRcutFromAtomType(AtomType* at) {
448      /**@todo */
449      GenericData* data;
450 <    double rcut = 0.0;
451 <
450 >    RealType rcut = 0.0;
451 >    
452      if (at->isLennardJones()) {
453 <        data = at->getPropertyByName("LennardJones");
454 <        if (data != NULL) {
455 <            LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
456 <
457 <            if (ljData != NULL) {
458 <                LJParam ljParam = ljData->getData();
459 <
460 <                //by default use 2.5*sigma as cutoff radius
461 <                rcut = 2.5 * ljParam.sigma;
462 <                
463 <            } else {
464 <                    sprintf( painCave.errMsg,
465 <                           "Can not cast GenericData to LJParam\n");
466 <                    painCave.severity = OOPSE_ERROR;
467 <                    painCave.isFatal = 1;
468 <                    simError();          
469 <            }            
470 <        } else {
471 <            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
472 <            painCave.severity = OOPSE_ERROR;
473 <            painCave.isFatal = 1;
474 <            simError();          
475 <        }
453 >      data = at->getPropertyByName("LennardJones");
454 >      if (data != NULL) {
455 >        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
456 >        
457 >        if (ljData != NULL) {
458 >          LJParam ljParam = ljData->getData();
459 >          
460 >          //by default use 2.5*sigma as cutoff radius
461 >          rcut = 2.5 * ljParam.sigma;
462 >          
463 >        } else {
464 >          sprintf( painCave.errMsg,
465 >                   "Can not cast GenericData to LJParam\n");
466 >          painCave.severity = OOPSE_ERROR;
467 >          painCave.isFatal = 1;
468 >          simError();          
469 >        }            
470 >      } else {
471 >        sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
472 >        painCave.severity = OOPSE_ERROR;
473 >        painCave.isFatal = 1;
474 >        simError();          
475 >      }
476      }
213 –
477      return rcut;    
478 < }
478 >  }
479 >  
480  
481 <
218 < ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
481 >  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
482      std::string forceFieldFilename(filename);
483      ifstrstream* ffStream = new ifstrstream();
484      
# Line 223 | Line 486 | ifstrstream* ForceField::openForceFieldFile(const std:
486      ffStream->open(forceFieldFilename.c_str());
487      if(!ffStream->is_open()){
488  
489 <        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
490 <        ffStream->open( forceFieldFilename.c_str() );
489 >      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
490 >      ffStream->open( forceFieldFilename.c_str() );
491  
492 <        //if current directory does not contain the force field file,
493 <        //try to open it in the path        
494 <        if(!ffStream->is_open()){
492 >      //if current directory does not contain the force field file,
493 >      //try to open it in the path        
494 >      if(!ffStream->is_open()){
495  
496 <            sprintf( painCave.errMsg,
497 <               "Error opening the force field parameter file:\n"
498 <               "\t%s\n"
499 <               "\tHave you tried setting the FORCE_PARAM_PATH environment "
500 <               "variable?\n",
501 <               forceFieldFilename.c_str() );
502 <            painCave.severity = OOPSE_ERROR;
503 <            painCave.isFatal = 1;
504 <            simError();
505 <        }
496 >        sprintf( painCave.errMsg,
497 >                 "Error opening the force field parameter file:\n"
498 >                 "\t%s\n"
499 >                 "\tHave you tried setting the FORCE_PARAM_PATH environment "
500 >                 "variable?\n",
501 >                 forceFieldFilename.c_str() );
502 >        painCave.severity = OOPSE_ERROR;
503 >        painCave.isFatal = 1;
504 >        simError();
505 >      }
506      }  
244 –
507      return ffStream;
508 +  }
509  
510 < }
511 <
510 >  void ForceField::setFortranForceOptions(){
511 >    ForceOptions theseFortranOptions;
512 >    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
513 >    setfForceOptions(&theseFortranOptions);
514 >  }
515   } //end namespace oopse

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines