--- trunk/src/primitives/Inversion.cpp	2009/11/25 20:02:06	1390
+++ trunk/src/primitives/Inversion.cpp	2013/06/18 21:06:25	1890
@@ -35,10 +35,14 @@
  *                                                                      
  * [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/Inversion.hpp"
 
 namespace OpenMD {
@@ -46,9 +50,11 @@ namespace OpenMD {
   Inversion::Inversion(Atom *atom1, Atom *atom2, Atom *atom3, 
 		       Atom *atom4, InversionType *it) :
     atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), 
-    inversionType_(it) { }
+    inversionType_(it) { 
+    inversionKey_ = inversionType_->getKey();
+  }
   
-  void Inversion::calcForce(RealType& angle) {
+  void Inversion::calcForce(RealType& angle, bool doParticlePot) {
     
     // In OpenMD's version of an inversion, the central atom
     // comes first.  However, to get the planarity in a typical cosine
@@ -82,7 +88,22 @@ namespace OpenMD {
     if (cos_phi < -1.0) cos_phi = -1.0;
 
     RealType dVdcosPhi;
-    inversionType_->calcForce(cos_phi, potential_, dVdcosPhi);
+    switch (inversionKey_) {
+    case itCosAngle:
+      inversionType_->calcForce(cos_phi, potential_, dVdcosPhi);
+      break;
+    case itAngle:
+      RealType phi = acos(cos_phi);
+      RealType dVdPhi;
+      inversionType_->calcForce(phi, potential_, dVdPhi);
+      RealType sin_phi = sqrt(1.0 - cos_phi * cos_phi);   
+      if (fabs(sin_phi) < 1.0E-6) {
+        sin_phi = 1.0E-6;
+      }
+      dVdcosPhi = dVdPhi / sin_phi;
+      break;
+    }
+    
     Vector3d f1 ;
     Vector3d f2 ;
     Vector3d f3 ;
@@ -109,11 +130,13 @@ namespace OpenMD {
     atom4_->addFrc(-f2);
     atom3_->addFrc(-f3);
 
-    atom1_->addParticlePot(potential_);
-    atom2_->addParticlePot(potential_);
-    atom3_->addParticlePot(potential_);
-    atom4_->addParticlePot(potential_);
-
+    if (doParticlePot) { 
+      atom1_->addParticlePot(potential_);
+      atom2_->addParticlePot(potential_);
+      atom3_->addParticlePot(potential_);
+      atom4_->addParticlePot(potential_);
+    }
+    
     angle = acos(cos_phi) /M_PI * 180.0;
   }