--- trunk/src/math/RectMatrix.hpp	2010/05/10 17:28:26	1442
+++ trunk/src/math/RectMatrix.hpp	2013/08/05 21:46:11	1924
@@ -35,8 +35,9 @@
  *                                                                      
  * [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).
  */
  
 /**
@@ -220,7 +221,7 @@ namespace OpenMD {
     /**
      * Tests if this matrix is identical to matrix m
      * @return true if this matrix is equal to the matrix m, return false otherwise
-     * @m matrix to be compared
+     * @param m matrix to be compared
      *
      * @todo replace operator == by template function equal
      */
@@ -236,7 +237,7 @@ namespace OpenMD {
     /**
      * Tests if this matrix is not equal to matrix m
      * @return true if this matrix is not equal to the matrix m, return false otherwise
-     * @m matrix to be compared
+     * @param m matrix to be compared
      */
     bool operator !=(const RectMatrix<Real, Row, Col>& m) {
       return !(*this == m);
@@ -506,7 +507,7 @@ namespace OpenMD {
   }
     
   /**
-   * Return the multiplication of  a matrix and a vector  (m * v). 
+   * Returns the multiplication of  a matrix and a vector  (m * v). 
    * @return the multiplication of a matrix and a vector
    * @param m the matrix
    * @param v the vector
@@ -518,6 +519,23 @@ namespace OpenMD {
     for (unsigned int i = 0; i < Row ; i++)
       for (unsigned int j = 0; j < Col ; j++)            
 	result[i] += m(i, j) * v[j];
+            
+    return result;                                                                 
+  }
+
+  /**
+   * Returns the multiplication of a vector transpose and a matrix  (v^T * m). 
+   * @return the multiplication of a vector transpose and a matrix
+   * @param v the vector
+   * @param m the matrix
+   */
+  template<typename Real, unsigned int Row, unsigned int Col>
+  inline Vector<Real, Col> operator *(const Vector<Real, Row>& v, const RectMatrix<Real, Row, Col>& m) {
+    Vector<Real, Row> result;
+    
+    for (unsigned int i = 0; i < Col ; i++)
+      for (unsigned int j = 0; j < Row ; j++)            
+	result[i] += v[j] * m(j, i);
             
     return result;                                                                 
   }
@@ -537,7 +555,72 @@ namespace OpenMD {
     return result;
   }    
 
+    
+    /**
+     * Returns the tensor contraction (double dot product) of two rank 2
+     * tensors (or Matrices)
+     *
+     * \f[ \mathbf{A} \colon \! \mathbf{B} = \sum_\alpha \sum_\beta \mathbf{A}_{\alpha \beta} B_{\alpha \beta} \f] 
+     *
+     * @param t1 first tensor
+     * @param t2 second tensor
+     * @return the tensor contraction (double dot product) of t1 and t2
+     */
+  template<typename Real, unsigned int Row, unsigned int Col> 
+  inline Real doubleDot( const RectMatrix<Real, Row, Col>& t1, 
+                         const RectMatrix<Real, Row, Col>& t2 ) {
+    Real tmp;
+    tmp = 0;
+    
+    for (unsigned int i = 0; i < Row; i++)
+      for (unsigned int j =0; j < Col; j++)
+        tmp += t1(i,j) * t2(i,j);
+    
+    return tmp;
+  }
+  
+
+  
   /**
+   * Returns the vector (cross) product of two matrices.  This
+   * operation is defined in:
+   *
+   * W. Smith, "Point Multipoles in the Ewald Summation (Revisited),"
+   * CCP5 Newsletter No 46., pp. 18-30.
+   *
+   * Equation 21 defines:
+   * \f[
+   * V_alpha = \sum_\beta \left[ A_{\alpha+1,\beta} * B_{\alpha+2,\beta} 
+                           -A_{\alpha+2,\beta} * B_{\alpha+2,\beta} \right]
+   * \f]
+
+   * where \f[\alpha+1\f] and \f[\alpha+2\f] are regarded as cyclic
+   * permuations of the matrix indices (i.e. for a 3x3 matrix, when
+   * \f[\alpha = 2\f], \f[\alpha + 1 = 3 \f], and \f[\alpha + 2 = 1 \f] ).
+   *
+   * @param t1 first matrix
+   * @param t2 second matrix
+   * @return the cross product (vector product) of t1 and t2
+   */
+  template<typename Real, unsigned int Row, unsigned int Col>
+  inline Vector<Real, Row> cross( const RectMatrix<Real, Row, Col>& t1, 
+                                  const RectMatrix<Real, Row, Col>& t2 ) {
+    Vector<Real, Row> result;
+    unsigned int i1;
+    unsigned int i2;
+    
+    for (unsigned int i = 0; i < Row; i++) {
+      i1 = (i+1)%Row;
+      i2 = (i+2)%Row;
+      for (unsigned int j = 0; j < Col; j++) {
+        result[i] += t1(i1,j) * t2(i2,j) - t1(i2,j) * t2(i1,j);
+      }
+    }    
+    return result;
+  }
+  
+  
+  /**
    * Write to an output stream
    */
   template<typename Real,  unsigned int Row, unsigned int Col>