--- trunk/src/nonbonded/Morse.cpp	2013/06/16 15:15:42	1879
+++ trunk/src/nonbonded/Morse.cpp	2014/09/02 18:31:44	2017
@@ -56,10 +56,16 @@ namespace OpenMD {
   
   void Morse::initialize() {    
 
+    Mtypes.clear();
+    Mtids.clear();
+    MixingMap.clear();
+    Mtids.resize( forceField_->getNAtomType(), -1);
+
     ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
     ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
-    NonBondedInteractionType* nbt;
     ForceField::NonBondedInteractionTypeContainer::KeyType keys;
+    NonBondedInteractionType* nbt;
+    int mtid1, mtid2;
 
     for (nbt = nbiTypes->beginType(j); nbt != NULL; 
          nbt = nbiTypes->nextType(j)) {
@@ -69,6 +75,19 @@ namespace OpenMD {
         AtomType* at1 = forceField_->getAtomType(keys[0]);
         AtomType* at2 = forceField_->getAtomType(keys[1]);
 
+        int atid1 = at1->getIdent();
+        if (Mtids[atid1] == -1) {
+          mtid1 = Mtypes.size();
+          Mtypes.insert(atid1);
+          Mtids[atid1] = mtid1;         
+        }
+        int atid2 = at2->getIdent();
+        if (Mtids[atid2] == -1) {
+          mtid2 = Mtypes.size();
+          Mtypes.insert(atid2);
+          Mtids[atid2] = mtid2;
+        }
+    
         MorseInteractionType* mit = dynamic_cast<MorseInteractionType*>(nbt);
 
         if (mit == NULL) {
@@ -103,13 +122,17 @@ namespace OpenMD {
     mixer.beta = beta;
     mixer.variant = mt;
 
-    pair<AtomType*, AtomType*> key1, key2;
-    key1 = make_pair(atype1, atype2);
-    key2 = make_pair(atype2, atype1);
+    int mtid1 = Mtids[atype1->getIdent()];
+    int mtid2 = Mtids[atype2->getIdent()];
+    int nM = Mtypes.size();
+
+    MixingMap.resize(nM);
+    MixingMap[mtid1].resize(nM);
     
-    MixingMap[key1] = mixer;
-    if (key2 != key1) {
-      MixingMap[key2] = mixer;
+    MixingMap[mtid1][mtid2] = mixer;
+    if (mtid2 != mtid1) {
+      MixingMap[mtid2].resize(nM);
+      MixingMap[mtid2][mtid1] = mixer;
     }    
   }
   
@@ -117,104 +140,104 @@ namespace OpenMD {
 
     if (!initialized_) initialize();
     
-    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
-    it = MixingMap.find( idat.atypes );
-    if (it != MixingMap.end()) {
-      MorseInteractionData mixer = (*it).second;
+    MorseInteractionData &mixer = MixingMap[Mtids[idat.atid1]][Mtids[idat.atid2]];
+
+    RealType myPot = 0.0;
+    RealType myPotC = 0.0;
+    RealType myDeriv = 0.0;
+    RealType myDerivC = 0.0;
+    
+    RealType De = mixer.De;
+    RealType Re = mixer.Re;
+    RealType beta = mixer.beta;
+    MorseType variant = mixer.variant;
+    
+    // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
+    
+    RealType expt     = -beta*( *(idat.rij) - Re);
+    RealType expfnc   = exp(expt);
+    RealType expfnc2  = expfnc*expfnc;
+    
+    RealType exptC = 0.0;
+    RealType expfncC = 0.0;
+    RealType expfnc2C = 0.0;
+    
+    if (idat.shiftedPot || idat.shiftedForce) {
+      exptC     = -beta*( *(idat.rcut) - Re);
+      expfncC   = exp(exptC);
+      expfnc2C  = expfncC*expfncC;
+    }
+    
+    
+    switch(variant) {
+    case mtShifted : {
       
-      RealType myPot = 0.0;
-      RealType myPotC = 0.0;
-      RealType myDeriv = 0.0;
-      RealType myDerivC = 0.0;
+      myPot  = De * (expfnc2  - 2.0 * expfnc);
+      myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
       
-      RealType De = mixer.De;
-      RealType Re = mixer.Re;
-      RealType beta = mixer.beta;
-      MorseType variant = mixer.variant;
-      
-      // V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2)
-      
-      RealType expt     = -beta*( *(idat.rij) - Re);
-      RealType expfnc   = exp(expt);
-      RealType expfnc2  = expfnc*expfnc;
-      
-      RealType exptC = 0.0;
-      RealType expfncC = 0.0;
-      RealType expfnc2C = 0.0;
-      
-      if (idat.shiftedPot || idat.shiftedForce) {
-        exptC     = -beta*( *(idat.rcut) - Re);
-        expfncC   = exp(exptC);
-        expfnc2C  = expfncC*expfncC;
+      if (idat.shiftedPot) {
+        myPotC = De * (expfnc2C - 2.0 * expfncC);
+        myDerivC = 0.0;
+      } else if (idat.shiftedForce) {
+        myPotC = De * (expfnc2C - 2.0 * expfncC);
+        myDerivC  = 2.0 * De * beta * (expfncC - expfnc2C);
+        myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
+      } else {
+        myPotC = 0.0;
+        myDerivC = 0.0;
       }
       
-      
-      switch(variant) {
-      case mtShifted : {
+      break;
+    }
+    case mtRepulsive : {
         
-        myPot  = De * (expfnc2  - 2.0 * expfnc);
-        myDeriv   = 2.0 * De * beta * (expfnc - expfnc2);
-        
-        if (idat.shiftedPot) {
-          myPotC = De * (expfnc2C - 2.0 * expfncC);
-          myDerivC = 0.0;
-        } else if (idat.shiftedForce) {
-          myPotC = De * (expfnc2C - 2.0 * expfncC);
-          myDerivC  = 2.0 * De * beta * (expfncC - expfnc2C);
-          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut) );
-        } else {
-          myPotC = 0.0;
-          myDerivC = 0.0;
-        }
-        
-        break;
+      myPot  = De * expfnc2;
+      myDeriv  = -2.0 * De * beta * expfnc2;
+      
+      if (idat.shiftedPot) {
+        myPotC = De * expfnc2C;
+        myDerivC = 0.0;
+      } else if (idat.shiftedForce) {
+        myPotC = De * expfnc2C;
+        myDerivC = -2.0 * De * beta * expfnc2C;
+        myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
+      } else {
+        myPotC = 0.0;
+        myDerivC = 0.0;
       }
-      case mtRepulsive : {
-        
-        myPot  = De * expfnc2;
-        myDeriv  = -2.0 * De * beta * expfnc2;
-        
-        if (idat.shiftedPot) {
-          myPotC = De * expfnc2C;
-          myDerivC = 0.0;
-        } else if (idat.shiftedForce) {
-          myPotC = De * expfnc2C;
-          myDerivC = -2.0 * De * beta * expfnc2C;
-          myPotC += myDerivC * ( *(idat.rij) - *(idat.rcut));
-        } else {
-          myPotC = 0.0;
-          myDerivC = 0.0;
-        }
-        
-        break;
-      }
-      case mtUnknown: {
-        // don't know what to do so don't do anything
-        break;
-      }
-      }
       
-      RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
-      *(idat.vpair) += pot_temp;
-      
-      RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
-      
-      (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
-      *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
+      break;
     }
-    return;
+    case mtUnknown: {
+      // don't know what to do so don't do anything
+      break;
+    }
+    }
     
-  }
+
+    RealType pot_temp = *(idat.vdwMult) * (myPot - myPotC);
+    *(idat.vpair) += pot_temp;
     
+    RealType dudr = *(idat.sw) * *(idat.vdwMult) * (myDeriv - myDerivC);
+
+    
+    (*(idat.pot))[VANDERWAALS_FAMILY] += *(idat.sw) * pot_temp;
+    *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
+    
+    return;    
+  }
+
   RealType Morse::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
     if (!initialized_) initialize();   
-    map<pair<AtomType*, AtomType*>, MorseInteractionData>::iterator it;
-    it = MixingMap.find(atypes);
-    if (it == MixingMap.end()) 
-      return 0.0;
-    else  {
-      MorseInteractionData mixer = (*it).second;
-
+    
+    int atid1 = atypes.first->getIdent();
+    int atid2 = atypes.second->getIdent();
+    int mtid1 = Mtids[atid1];
+    int mtid2 = Mtids[atid2];
+    
+    if ( mtid1 == -1 || mtid2 == -1) return 0.0;
+    else {      
+      MorseInteractionData mixer = MixingMap[mtid1][mtid2];
       RealType Re = mixer.Re;
       RealType beta = mixer.beta;
       // This value of the r corresponds to an energy about 1.48% of