--- trunk/src/primitives/Torsion.cpp	2009/11/25 20:02:06	1390
+++ trunk/src/primitives/Torsion.cpp	2013/12/05 18:19:26	1953
@@ -35,24 +35,34 @@
  *                                                                      
  * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).             
  * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
- * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
- * [4]  Vardeman & Gezelter, in progress (2009).                        
+ * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
+ * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
+ * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
  */
  
+#include "config.h"
+#include <cmath>
+
 #include "primitives/Torsion.hpp"
 
 namespace OpenMD {
 
   Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
-		   TorsionType *tt) :
-    atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
+		   TorsionType *tt) : ShortRangeInteraction(),
+                                      torsionType_(tt) {
+    atoms_.resize(4);
+    atoms_[0] = atom1;
+    atoms_[1] = atom2;
+    atoms_[2] = atom3;
+    atoms_[3] = atom4;
+  }
 
-  void Torsion::calcForce(RealType& angle) {
+  void Torsion::calcForce(RealType& angle, bool doParticlePot) {
 
-    Vector3d pos1 = atom1_->getPos();
-    Vector3d pos2 = atom2_->getPos();
-    Vector3d pos3 = atom3_->getPos();
-    Vector3d pos4 = atom4_->getPos();
+    Vector3d pos1 = atoms_[0]->getPos();
+    Vector3d pos2 = atoms_[1]->getPos();
+    Vector3d pos3 = atoms_[2]->getPos();
+    Vector3d pos4 = atoms_[3]->getPos();
 
     Vector3d r21 = pos1 - pos2;
     Vector3d r32 = pos2 - pos3;
@@ -63,42 +73,48 @@ namespace OpenMD {
     RealType rA = A.length();
     Vector3d B = cross(r32, r43);
     RealType rB = B.length();
-    Vector3d C = cross(r32, A);
-    RealType rC = C.length();
 
+    /* 
+       If either of the two cross product vectors is tiny, that means
+       the three atoms involved are colinear, and the torsion angle is
+       going to be undefined.  The easiest check for this problem is
+       to use the product of the two lengths.
+    */
+    if (rA * rB < OpenMD::epsilon) return;
+    
     A.normalize();
-    B.normalize();
-    C.normalize();
+    B.normalize();  
     
     //  Calculate the sin and cos
     RealType cos_phi = dot(A, B) ;
     if (cos_phi > 1.0) cos_phi = 1.0;
     if (cos_phi < -1.0) cos_phi = -1.0; 
-
+    
     RealType dVdcosPhi;
     torsionType_->calcForce(cos_phi, potential_, dVdcosPhi);
     Vector3d f1 ;
     Vector3d f2 ;
     Vector3d f3 ;
-
+    
     Vector3d dcosdA = (cos_phi * A - B) /rA;
     Vector3d dcosdB = (cos_phi * B - A) /rB;
-
+    
     f1 = dVdcosPhi * cross(r32, dcosdA);
     f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
     f3 = dVdcosPhi * cross(dcosdB, r32);
     
-    atom1_->addFrc(f1);
-    atom2_->addFrc(f2 - f1);
-    atom3_->addFrc(f3 - f2);
-    atom4_->addFrc(-f3);
-
-    atom1_->addParticlePot(potential_);
-    atom2_->addParticlePot(potential_);
-    atom3_->addParticlePot(potential_);
-    atom4_->addParticlePot(potential_);
-
-    angle = acos(cos_phi) /M_PI * 180.0;
-  }
-
+    atoms_[0]->addFrc(f1);
+    atoms_[1]->addFrc(f2 - f1);
+    atoms_[2]->addFrc(f3 - f2);
+    atoms_[3]->addFrc(-f3);
+    
+    if (doParticlePot) {
+      atoms_[0]->addParticlePot(potential_);
+      atoms_[1]->addParticlePot(potential_);
+      atoms_[2]->addParticlePot(potential_);
+      atoms_[3]->addParticlePot(potential_);
+    }
+    
+    angle = acos(cos_phi) /M_PI * 180.0;    
+  }  
 }