[ViewVC] Diff of: OpenMD/trunk/src/brains/ForceField.cpp

Comparing trunk/src/UseTheForce/ForceField.cpp (file contents):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 1282 by gezelter, Wed Jul 30 18:11:19 2008 UTC

+ * @version 1.0
+ */
-+
+#include <algorithm>
+#include "UseTheForce/ForceField.hpp"
+#include "utils/simError.h"
-+
+#include "utils/Tuple.hpp"
+#include "UseTheForce/DarkSide/atype_interface.h"
+#include "UseTheForce/DarkSide/fForceOptions_interface.h"
+#include "UseTheForce/DarkSide/switcheroo_interface.h"
+    return atomTypeCont_.find(keys);
-<
+  BondType* ForceField::getBondType(const std::string &at1, const std::string &at2) {
->
+  BondType* ForceField::getBondType(const std::string &at1,
->
+                                    const std::string &at2) {
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
+    if (bondType) {
+      return bondType;
+    } else {
-<
+      //if no exact match found, try wild card match
-<
+      return bondTypeCont_.find(keys, wildCardAtomTypeName_);
-<
+    }
->
+      AtomType* atype1;
->
+      AtomType* atype2;
->
+      std::vector<std::string> at1key;
->
+      at1key.push_back(at1);
->
+      atype1 = atomTypeCont_.find(at1key);
->
->
+      std::vector<std::string> at2key;
->
+      at2key.push_back(at2);
->
+      atype2 = atomTypeCont_.find(at2key);
-<
+  }
->
+      // query atom types for their chains of responsibility
->
+      std::vector<AtomType*> at1Chain = atype1->allYourBase();
->
+      std::vector<AtomType*> at2Chain = atype2->allYourBase();
-<
+  BendType* ForceField::getBendType(const std::string &at1, const std::string &at2,
->
+      std::vector<AtomType*>::iterator i;
->
+      std::vector<AtomType*>::iterator j;
->
->
+      int ii = 0;
->
+      int jj = 0;
->
+      int bondTypeScore;
->
->
+      std::vector<std::pair<int, std::vector<std::string> > > foundBonds;
->
->
+      for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
->
+        jj = 0;
->
+        for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
->
->
+          bondTypeScore = ii + jj;
->
->
+          std::vector<std::string> myKeys;
->
+          myKeys.push_back((*i)->getName());
->
+          myKeys.push_back((*j)->getName());
->
->
+          std::cerr << "looking for " << myKeys[0] << " " << myKeys[1] << "\n";
->
+          BondType* bondType = bondTypeCont_.find(myKeys);
->
+          if (bondType) {
->
+            foundBonds.push_back(std::make_pair(bondTypeScore, myKeys));
->
+          }
->
+          jj++;
->
+        }
->
+        ii++;
->
+      }
->
->
->
+      if (foundBonds.size() > 0) {
->
+        // sort the foundBonds by the score:
->
+        std::sort(foundBonds.begin(), foundBonds.end());
->
->
+        int bestScore = foundBonds[0].first;
->
+        std::vector<std::string> theKeys = foundBonds[0].second;
->
->
+        std::cout << "best matching bond = " << theKeys[0] << "\t" << theKeys[1]  << "\t(score = "<< bestScore << ")\n";
->
+        BondType* bestType = bondTypeCont_.find(theKeys);
->
->
+        return bestType;
->
+      } else {
->
+        //if no exact match found, try wild card match
->
+        return bondTypeCont_.find(keys, wildCardAtomTypeName_);
->
+      }
->
+    }
->
+  }
->
->
+  BendType* ForceField::getBendType(const std::string &at1,
->
+                                    const std::string &at2,
+                                    const std::string &at3) {
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    if (bendType) {
+      return bendType;
+    } else {
-<
+      //if no exact match found, try wild card match
-<
+      return bendTypeCont_.find(keys, wildCardAtomTypeName_);
->
->
+      AtomType* atype1;
->
+      AtomType* atype2;
->
+      AtomType* atype3;
->
+      std::vector<std::string> at1key;
->
+      at1key.push_back(at1);
->
+      atype1 = atomTypeCont_.find(at1key);
->
->
+      std::vector<std::string> at2key;
->
+      at2key.push_back(at2);
->
+      atype2 = atomTypeCont_.find(at2key);
->
->
+      std::vector<std::string> at3key;
->
+      at3key.push_back(at3);
->
+      atype3 = atomTypeCont_.find(at3key);
->
->
+      // query atom types for their chains of responsibility
->
+      std::vector<AtomType*> at1Chain = atype1->allYourBase();
->
+      std::vector<AtomType*> at2Chain = atype2->allYourBase();
->
+      std::vector<AtomType*> at3Chain = atype3->allYourBase();
->
->
+      std::vector<AtomType*>::iterator i;
->
+      std::vector<AtomType*>::iterator j;
->
+      std::vector<AtomType*>::iterator k;
->
->
+      int ii = 0;
->
+      int jj = 0;
->
+      int kk = 0;
->
+      int IKscore;
->
->
+      std::vector<tuple3<int, int, std::vector<std::string> > > foundBends;
->
->
+      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
->
+        ii = 0;
->
+        for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
->
+          kk = 0;
->
+          for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
->
->
+            IKscore = ii + kk;
->
->
+            std::vector<std::string> myKeys;
->
+            myKeys.push_back((*i)->getName());
->
+            myKeys.push_back((*j)->getName());
->
+            myKeys.push_back((*k)->getName());
->
->
+            BendType* bendType = bendTypeCont_.find(myKeys);
->
+            if (bendType) {
->
+              foundBends.push_back( make_tuple3(jj, IKscore, myKeys) );
->
+            }
->
+            kk++;
->
+          }
->
+          ii++;
->
+        }
->
+        jj++;
->
+      }
->
->
+      if (foundBends.size() > 0) {
->
+        std::sort(foundBends.begin(), foundBends.end());
->
+        int jscore = foundBends[0].first;
->
+        int ikscore = foundBends[0].second;
->
+        std::vector<std::string> theKeys = foundBends[0].third;
->
->
+        std::cout << "best matching bend = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t(scores = "<< jscore << "\t" << ikscore << ")\n";
->
->
+        BendType* bestType = bendTypeCont_.find(theKeys);
->
+        return bestType;
->
+      } else {
->
+        //if no exact match found, try wild card match
->
+        return bendTypeCont_.find(keys, wildCardAtomTypeName_);
->
+      }
-<
+  TorsionType* ForceField::getTorsionType(const std::string &at1, const std::string &at2,
-<
+                                          const std::string &at3, const std::string &at4) {
->
+  TorsionType* ForceField::getTorsionType(const std::string &at1,
->
+                                          const std::string &at2,
->
+                                          const std::string &at3,
->
+                                          const std::string &at4) {
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
+    keys.push_back(at3);
+    keys.push_back(at4);
-+
-+
+    //try exact match first
+    TorsionType* torsionType = torsionTypeCont_.find(keys);
+    if (torsionType) {
+      return torsionType;
+    } else {
-–
+      //if no exact match found, try wild card match
-–
+      return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
-–
+    }
-–
-–
+    return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
-+
+      AtomType* atype1;
-+
+      AtomType* atype2;
-+
+      AtomType* atype3;
-+
+      AtomType* atype4;
-+
+      std::vector<std::string> at1key;
-+
+      at1key.push_back(at1);
-+
+      atype1 = atomTypeCont_.find(at1key);
-+
-+
+      std::vector<std::string> at2key;
-+
+      at2key.push_back(at2);
-+
+      atype2 = atomTypeCont_.find(at2key);
-+
-+
+      std::vector<std::string> at3key;
-+
+      at3key.push_back(at3);
-+
+      atype3 = atomTypeCont_.find(at3key);
-+
-+
+      std::vector<std::string> at4key;
-+
+      at4key.push_back(at4);
-+
+      atype4 = atomTypeCont_.find(at4key);
-+
-+
+      // query atom types for their chains of responsibility
-+
+      std::vector<AtomType*> at1Chain = atype1->allYourBase();
-+
+      std::vector<AtomType*> at2Chain = atype2->allYourBase();
-+
+      std::vector<AtomType*> at3Chain = atype3->allYourBase();
-+
+      std::vector<AtomType*> at4Chain = atype4->allYourBase();
-+
-+
+      std::vector<AtomType*>::iterator i;
-+
+      std::vector<AtomType*>::iterator j;
-+
+      std::vector<AtomType*>::iterator k;
-+
+      std::vector<AtomType*>::iterator l;
-+
-+
+      int ii = 0;
-+
+      int jj = 0;
-+
+      int kk = 0;
-+
+      int ll = 0;
-+
+      int ILscore;
-+
+      int JKscore;
-+
-+
+      std::vector<tuple3<int, int, std::vector<std::string> > > foundTorsions;
-+
-+
+      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
-+
+        kk = 0;
-+
+        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
-+
+          ii = 0;
-+
+          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
-+
+            ll = 0;
-+
+            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
-+
-+
+              ILscore = ii + ll;
-+
+              JKscore = jj + kk;
-+
-+
+              std::vector<std::string> myKeys;
-+
+              myKeys.push_back((*i)->getName());
-+
+              myKeys.push_back((*j)->getName());
-+
+              myKeys.push_back((*k)->getName());
-+
+              myKeys.push_back((*l)->getName());
-+
-+
+              TorsionType* torsionType = torsionTypeCont_.find(myKeys);
-+
+              if (torsionType) {
-+
+                foundTorsions.push_back( make_tuple3(JKscore, ILscore, myKeys) );
-+
+              }
-+
+              ll++;
-+
+            }
-+
+            ii++;
-+
+          }
-+
+          kk++;
-+
+        }
-+
+        jj++;
-+
+      }
-+
-+
+      if (foundTorsions.size() > 0) {
-+
+        std::sort(foundTorsions.begin(), foundTorsions.end());
-+
+        int jkscore = foundTorsions[0].first;
-+
+        int ilscore = foundTorsions[0].second;
-+
+        std::vector<std::string> theKeys = foundTorsions[0].third;
-+
-+
+        std::cout << "best matching torsion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< jkscore << "\t" << ilscore << ")\n";
-+
-+
+        TorsionType* bestType = torsionTypeCont_.find(theKeys);
-+
+        return bestType;
-+
+      } else {
-+
+        //if no exact match found, try wild card match
-+
+        return torsionTypeCont_.find(keys, wildCardAtomTypeName_);
-+
+      }
-+
+    }
-<
+  BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){
->
+  InversionType* ForceField::getInversionType(const std::string &at1,
->
+                                              const std::string &at2,
->
+                                              const std::string &at3,
->
+                                              const std::string &at4) {
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
-+
+    keys.push_back(at3);
-+
+    keys.push_back(at4);
-+
-+
+    //try exact match first
-+
+    InversionType* inversionType = inversionTypeCont_.find(keys);
-+
+    if (inversionType) {
-+
+      return inversionType;
-+
+    } else {
-+
-+
+      AtomType* atype1;
-+
+      AtomType* atype2;
-+
+      AtomType* atype3;
-+
+      AtomType* atype4;
-+
+      std::vector<std::string> at1key;
-+
+      at1key.push_back(at1);
-+
+      atype1 = atomTypeCont_.find(at1key);
-+
-+
+      std::vector<std::string> at2key;
-+
+      at2key.push_back(at2);
-+
+      atype2 = atomTypeCont_.find(at2key);
-+
-+
+      std::vector<std::string> at3key;
-+
+      at3key.push_back(at3);
-+
+      atype3 = atomTypeCont_.find(at3key);
-+
-+
+      std::vector<std::string> at4key;
-+
+      at4key.push_back(at4);
-+
+      atype4 = atomTypeCont_.find(at4key);
-+
-+
+      // query atom types for their chains of responsibility
-+
+      std::vector<AtomType*> at1Chain = atype1->allYourBase();
-+
+      std::vector<AtomType*> at2Chain = atype2->allYourBase();
-+
+      std::vector<AtomType*> at3Chain = atype3->allYourBase();
-+
+      std::vector<AtomType*> at4Chain = atype4->allYourBase();
-+
-+
+      std::vector<AtomType*>::iterator i;
-+
+      std::vector<AtomType*>::iterator j;
-+
+      std::vector<AtomType*>::iterator k;
-+
+      std::vector<AtomType*>::iterator l;
-+
-+
+      int ii = 0;
-+
+      int jj = 0;
-+
+      int kk = 0;
-+
+      int ll = 0;
-+
+      int Iscore;
-+
+      int JKLscore;
-+
-+
+      std::vector<tuple3<int, int, std::vector<std::string> > > foundInversions;
-+
-+
+      for (j = at2Chain.begin(); j != at2Chain.end(); j++) {
-+
+        kk = 0;
-+
+        for (k = at3Chain.begin(); k != at3Chain.end(); k++) {
-+
+          ii = 0;
-+
+          for (i = at1Chain.begin(); i != at1Chain.end(); i++) {
-+
+            ll = 0;
-+
+            for (l = at4Chain.begin(); l != at4Chain.end(); l++) {
-+
-+
+              Iscore = ii;
-+
+              JKLscore = jj + kk + ll;
-+
-+
+              std::vector<std::string> myKeys;
-+
+              myKeys.push_back((*i)->getName());
-+
+              myKeys.push_back((*j)->getName());
-+
+              myKeys.push_back((*k)->getName());
-+
+              myKeys.push_back((*l)->getName());
-+
-+
+              InversionType* inversionType = inversionTypeCont_.find(myKeys);
-+
+              if (inversionType) {
-+
+                foundInversions.push_back( make_tuple3(Iscore, JKLscore, myKeys) );
-+
+              }
-+
+              ll++;
-+
+            }
-+
+            ii++;
-+
+          }
-+
+          kk++;
-+
+        }
-+
+        jj++;
-+
+      }
-+
-+
+      if (foundInversions.size() > 0) {
-+
+        std::sort(foundInversions.begin(), foundInversions.end());
-+
+        int iscore = foundInversions[0].first;
-+
+        int jklscore = foundInversions[0].second;
-+
+        std::vector<std::string> theKeys = foundInversions[0].third;
-+
-+
+        std::cout << "best matching inversion = " << theKeys[0] << "\t" <<theKeys[1]  << "\t" << theKeys[2] << "\t" << theKeys[3] << "\t(scores = "<< iscore << "\t" << jklscore << ")\n";
-+
-+
+        InversionType* bestType = inversionTypeCont_.find(theKeys);
-+
+        return bestType;
-+
+      } else {
-+
+        //if no exact match found, try wild card match
-+
+        return inversionTypeCont_.find(keys, wildCardAtomTypeName_);
-+
+      }
-+
+    }
-+
+  }
-+
-+
+  NonBondedInteractionType* ForceField::getNonBondedInteractionType(const std::string &at1, const std::string &at2) {
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
-+
-+
+    //try exact match first
-+
+    NonBondedInteractionType* nbiType = nonBondedInteractionTypeCont_.find(keys);
-+
+    if (nbiType) {
-+
+      return nbiType;
-+
+    } else {
-+
+      //if no exact match found, try wild card match
-+
+      return nonBondedInteractionTypeCont_.find(keys, wildCardAtomTypeName_);
-+
+    }
-+
+  }
-+
-+
+  BondType* ForceField::getExactBondType(const std::string &at1,
-+
+                                         const std::string &at2){
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
+    return bondTypeCont_.find(keys);
-<
-<
+  BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
->
->
+  BendType* ForceField::getExactBendType(const std::string &at1,
->
+                                         const std::string &at2,
+                                         const std::string &at3){
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at3);
+    return bendTypeCont_.find(keys);
-<
-<
+  TorsionType* ForceField::getExactTorsionType(const std::string &at1, const std::string &at2,
-<
+                                               const std::string &at3, const std::string &at4){
->
->
+  TorsionType* ForceField::getExactTorsionType(const std::string &at1,
->
+                                               const std::string &at2,
->
+                                               const std::string &at3,
->
+                                               const std::string &at4){
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
+    keys.push_back(at4);
+    return torsionTypeCont_.find(keys);
-+
-+
+  InversionType* ForceField::getExactInversionType(const std::string &at1,
-+
+                                                   const std::string &at2,
-+
+                                                   const std::string &at3,
-+
+                                                   const std::string &at4){
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
-+
+    keys.push_back(at3);
-+
+    keys.push_back(at4);
-+
+    return inversionTypeCont_.find(keys);
-+
+  }
-+
-+
+  NonBondedInteractionType* ForceField::getExactNonBondedInteractionType(const std::string &at1, const std::string &at2){
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
-+
+    return nonBondedInteractionTypeCont_.find(keys);
-+
+  }
-+
-+
+  bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
+    std::vector<std::string> keys;
+    keys.push_back(at);
+    return atomTypeCont_.add(keys, atomType);
-<
+  bool ForceField::addBondType(const std::string &at1, const std::string &at2, BondType* bondType) {
->
+  bool ForceField::replaceAtomType(const std::string &at, AtomType* atomType) {
+    std::vector<std::string> keys;
-+
+    keys.push_back(at);
-+
+    return atomTypeCont_.replace(keys, atomType);
-+
+  }
-+
-+
+  bool ForceField::addBondType(const std::string &at1, const std::string &at2,
-+
+                               BondType* bondType) {
-+
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
-<
+    return bondTypeCont_.add(keys, bondType);
-<
->
+    return bondTypeCont_.add(keys, bondType);
-<
->
+  bool ForceField::addBendType(const std::string &at1, const std::string &at2,
+                               const std::string &at3, BendType* bendType) {
+    std::vector<std::string> keys;
+    keys.push_back(at3);
+    return bendTypeCont_.add(keys, bendType);
-<
-<
+  bool ForceField::addTorsionType(const std::string &at1, const std::string &at2,
-<
+                                  const std::string &at3, const std::string &at4, TorsionType* torsionType) {
->
->
+  bool ForceField::addTorsionType(const std::string &at1,
->
+                                  const std::string &at2,
->
+                                  const std::string &at3,
->
+                                  const std::string &at4,
->
+                                  TorsionType* torsionType) {
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);
+    return torsionTypeCont_.add(keys, torsionType);
-+
+  bool ForceField::addInversionType(const std::string &at1,
-+
+                                    const std::string &at2,
-+
+                                    const std::string &at3,
-+
+                                    const std::string &at4,
-+
+                                    InversionType* inversionType) {
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
-+
+    keys.push_back(at3);
-+
+    keys.push_back(at4);
-+
+    return inversionTypeCont_.add(keys, inversionType);
-+
+  }
-+
-+
+  bool ForceField::addNonBondedInteractionType(const std::string &at1,
-+
+                                               const std::string &at2,
-+
+                                               NonBondedInteractionType* nbiType) {
-+
+    std::vector<std::string> keys;
-+
+    keys.push_back(at1);
-+
+    keys.push_back(at2);
-+
+    return nonBondedInteractionTypeCont_.add(keys, nbiType);
-+
+  }
-+
+  RealType ForceField::getRcutFromAtomType(AtomType* at) {
+    /**@todo */
+    GenericData* data;
+    RealType rcut = 0.0;
-<
->
+    if (at->isLennardJones()) {
+      data = at->getPropertyByName("LennardJones");
+      if (data != NULL) {
+        LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
-<
->
+        if (ljData != NULL) {
+          LJParam ljParam = ljData->getData();
-<
->
+          //by default use 2.5*sigma as cutoff radius
+          rcut = 2.5 * ljParam.sigma;
-<
->
+        } else {
+          sprintf( painCave.errMsg,
+                   "Can not cast GenericData to LJParam\n");
+        simError();
-–
+    return rcut;
-+
-–
+  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
+    std::string forceFieldFilename(filename);
+    ifstrstream* ffStream = new ifstrstream();
+        simError();
-–
+    return ffStream;
-–
+  void ForceField::setFortranForceOptions(){

Diff Legend

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

Comparing trunk/src/UseTheForce/ForceField.cpp (file contents): Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs. Revision 1282 by gezelter, Wed Jul 30 18:11:19 2008 UTC

Diff Legend

Comparing trunk/src/UseTheForce/ForceField.cpp (file contents):
Revision 963 by tim, Wed May 17 21:51:42 2006 UTC vs.
Revision 1282 by gezelter, Wed Jul 30 18:11:19 2008 UTC