--- trunk/src/UseTheForce/ForceField.cpp	2005/01/12 22:41:40	246
+++ trunk/src/UseTheForce/ForceField.cpp	2008/07/30 18:11:19	1282
@@ -1,4 +1,4 @@
- /*
+/*
  * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
  *
  * The University of Notre Dame grants you ("Licensee") a
@@ -39,37 +39,49 @@
  * such damages.
  */
  
- /**
-  * @file ForceField.cpp
-  * @author tlin
-  * @date 11/04/2004
-  * @time 22:51am
-  * @version 1.0
-  */
+/**
+ * @file ForceField.cpp
+ * @author tlin
+ * @date 11/04/2004
+ * @time 22:51am
+ * @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"
 namespace oopse {
 
-ForceField::ForceField() { 
+  ForceField::ForceField() { 
     char* tempPath; 
     tempPath = getenv("FORCE_PARAM_PATH");
 
     if (tempPath == NULL) {
-        //convert a macro from compiler to a string in c++
-        STR_DEFINE(ffPath_, FRC_PATH );
+      //convert a macro from compiler to a string in c++
+      STR_DEFINE(ffPath_, FRC_PATH );
     } else {
-        ffPath_ = tempPath;
+      ffPath_ = tempPath;
     }
-}
+  }
 
-AtomType* ForceField::getAtomType(const std::string &at) {
+
+  ForceField::~ForceField() {
+    deleteAtypes();
+    deleteSwitch();
+  }
+
+  AtomType* ForceField::getAtomType(const std::string &at) {
     std::vector<std::string> keys;
     keys.push_back(at);
     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);    
@@ -77,16 +89,73 @@ BondType* ForceField::getBondType(const std::string &a
     //try exact match first
     BondType* bondType = bondTypeCont_.find(keys);
     if (bondType) {
-        return 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,
-                        const std::string &at3) {
+      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);
     keys.push_back(at2);    
@@ -95,127 +164,452 @@ BendType* ForceField::getBendType(const std::string &a
     //try exact match first
     BendType* bendType = bendTypeCont_.find(keys);
     if (bendType) {
-        return 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;
+      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);
 
-BondType* ForceField::getExactBondType(const std::string &at1, const std::string &at2){ 
+      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_);
+      }
+    }
+  }
+
+  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);    
-    return bondTypeCont_.find(keys);
-}
+    keys.push_back(at3);    
+    keys.push_back(at4);    
 
-BendType* ForceField::getExactBendType(const std::string &at1, const std::string &at2,
-                            const std::string &at3){ 
+    //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,
+					 const std::string &at3){ 
+    std::vector<std::string> keys;
+    keys.push_back(at1);
+    keys.push_back(at2);    
     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(at3);    
     keys.push_back(at4);   
     return torsionTypeCont_.find(keys);
-}
-bool ForceField::addAtomType(const std::string &at, AtomType* atomType) {
+  }
+  
+  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);
-
-}
-
-bool ForceField::addBendType(const std::string &at1, const std::string &at2,
-                        const std::string &at3, BendType* bendType) {
+    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(at1);
     keys.push_back(at2);    
     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);    
     keys.push_back(at3);    
     keys.push_back(at4);    
     return torsionTypeCont_.add(keys, torsionType);
-}
+  }
 
-double ForceField::getRcutFromAtomType(AtomType* at) {
+  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;
-    double rcut = 0.0;
-
+    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");
-                    painCave.severity = OOPSE_ERROR;
-                    painCave.isFatal = 1;
-                    simError();          
-            }            
-        } else {
-            sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
-            painCave.severity = OOPSE_ERROR;
-            painCave.isFatal = 1;
-            simError();          
-        }
+      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");
+	  painCave.severity = OOPSE_ERROR;
+	  painCave.isFatal = 1;
+	  simError();          
+	}            
+      } else {
+	sprintf( painCave.errMsg, "Can not find Parameters for LennardJones\n");
+	painCave.severity = OOPSE_ERROR;
+	painCave.isFatal = 1;
+	simError();          
+      }
     }
-
     return rcut;    
-}
+  }
+  
 
-
-ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
+  ifstrstream* ForceField::openForceFieldFile(const std::string& filename) {
     std::string forceFieldFilename(filename);
     ifstrstream* ffStream = new ifstrstream();
     
@@ -223,27 +617,30 @@ ifstrstream* ForceField::openForceFieldFile(const std:
     ffStream->open(forceFieldFilename.c_str());
     if(!ffStream->is_open()){
 
-        forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
-        ffStream->open( forceFieldFilename.c_str() );
+      forceFieldFilename = ffPath_ + "/" + forceFieldFilename;
+      ffStream->open( forceFieldFilename.c_str() );
 
-        //if current directory does not contain the force field file,
-        //try to open it in the path        
-        if(!ffStream->is_open()){
+      //if current directory does not contain the force field file,
+      //try to open it in the path        
+      if(!ffStream->is_open()){
 
-            sprintf( painCave.errMsg,
-               "Error opening the force field parameter file:\n"
-               "\t%s\n"
-               "\tHave you tried setting the FORCE_PARAM_PATH environment "
-               "variable?\n",
-               forceFieldFilename.c_str() );
-            painCave.severity = OOPSE_ERROR;
-            painCave.isFatal = 1;
-            simError();
-        }
+	sprintf( painCave.errMsg,
+		 "Error opening the force field parameter file:\n"
+		 "\t%s\n"
+		 "\tHave you tried setting the FORCE_PARAM_PATH environment "
+		 "variable?\n",
+		 forceFieldFilename.c_str() );
+	painCave.severity = OOPSE_ERROR;
+	painCave.isFatal = 1;
+	simError();
+      }
     }  
-
     return ffStream;
+  }
 
-}
-
+  void ForceField::setFortranForceOptions(){
+    ForceOptions theseFortranOptions;
+    forceFieldOptions_.makeFortranOptions(theseFortranOptions);
+    setfForceOptions(&theseFortranOptions);
+  }
 } //end namespace oopse