? .cvsignore
Index: mortsrc/capbox/region.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/capbox/region.cpp,v
retrieving revision 1.1
diff -u -r1.1 region.cpp
--- mortsrc/capbox/region.cpp	15 Feb 2009 16:49:24 -0000	1.1
+++ mortsrc/capbox/region.cpp	27 Sep 2009 16:04:16 -0000
@@ -160,7 +160,7 @@
     {
         numvec reg = region( solute.getcord(), solute.getvdwr() );
         m_center = makevec( reg[0], reg[1], reg[2] );
-        m_extent = makevec( reg[3], reg[4], reg[5] ) + 2*shlext;
+        m_extent = makevec( reg[3], reg[4], reg[5] ) + scalar_numvec(3, 2*shlext);
 
         m_maxrad = *std::max_element( solute.getvdwr().begin(), solute.getvdwr().end() );
         m_clsnss = clsnss;
@@ -204,10 +204,10 @@
     {
         funstack_t::push( "and_region::and_region" );
 
-        numvec cnt1 = lhs.center();
+        numvec cnt1 = subvec(lhs.center(), 0, 3);
         numvec ext1 = lhs.extent();
 
-        numvec cnt2 = rhs.center();
+        numvec cnt2 = subvec(rhs.center(), 0, 3);
         numvec ext2 = rhs.extent();
 
         numvec min1 = cnt1 - 0.5*ext1;
Index: mortsrc/capbox/region.hpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/capbox/region.hpp,v
retrieving revision 1.1
diff -u -r1.1 region.hpp
--- mortsrc/capbox/region.hpp	15 Feb 2009 16:49:24 -0000	1.1
+++ mortsrc/capbox/region.hpp	27 Sep 2009 16:04:16 -0000
@@ -52,7 +52,7 @@
 
         // if no center is given, consider it is on origin
         box_region( const numvec& boxext )
-            : m_cnt(0.0, 3), m_ext(boxext), m_min(3), m_max(3)
+            : m_cnt(zero_numvec(3)), m_ext(boxext), m_min(3), m_max(3)
         {
             m_min = m_cnt - 0.5*boxext;
             m_max = m_cnt + 0.5*boxext;
@@ -83,7 +83,7 @@
         cap_region( const numvec& capcnt, double caprad )
            : m_cnt(capcnt),
              m_rad(caprad),
-             m_ext(2*caprad, 3)
+             m_ext(scalar_numvec(3, 2*caprad))
         {
         }
 
@@ -114,7 +114,7 @@
         }
 
         oct_region( const numvec& octext )
-            : m_cnt(0.0, 3), m_ext(octext)
+            : m_cnt(zero_numvec(3)), m_ext(octext)
         {
         }
 
@@ -153,9 +153,9 @@
 
         virtual int check_rndbox( const numvec& boxcnt, const numvec& boxext, double sphrad ) const;
 
-        virtual numvec center() const { return numvec(0.0, 3); }
+        virtual numvec center() const { return zero_numvec(3); }
 
-        virtual numvec extent() const { return numvec(DOUBLEMAX, 3); }
+        virtual numvec extent() const { return scalar_numvec(3, DOUBLEMAX); }
 
         int check_sphere( const numvec& sphcnt, double sphrad, double& mind2 ) const;
 
Index: mortsrc/capbox/solvate.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/capbox/solvate.cpp,v
retrieving revision 1.26
diff -u -r1.26 solvate.cpp
--- mortsrc/capbox/solvate.cpp	16 Feb 2009 03:31:39 -0000	1.26
+++ mortsrc/capbox/solvate.cpp	27 Sep 2009 16:04:16 -0000
@@ -131,7 +131,7 @@
         funstack_t::push( "solvatebox" );
 
         numvec extent = regionlize( solute, true );
-        extent += 2*buffer;
+        extent += scalar_numvec(extent.size(), 2*buffer);
 
         box_region box( extent );
         out_solute inn( solute, closeness );
Index: mortsrc/capbox/solvent.hpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/capbox/solvent.hpp,v
retrieving revision 1.1
diff -u -r1.1 solvent.hpp
--- mortsrc/capbox/solvent.hpp	15 Feb 2009 16:49:24 -0000	1.1
+++ mortsrc/capbox/solvent.hpp	27 Sep 2009 16:04:16 -0000
@@ -37,7 +37,7 @@
 
         const numvec center() const { return mort::center(*m_pmol); }
 
-        const numvec extent() const { return m_pmol->get_v(BOX)[std::slice(0,3,1)]; }
+        const numvec extent() const { return subvec(m_pmol->get_v(BOX), 0, 3); }
 
     private:
 
Index: mortsrc/common/crdalgo.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/common/crdalgo.cpp,v
retrieving revision 1.3
diff -u -r1.3 crdalgo.cpp
--- mortsrc/common/crdalgo.cpp	23 Feb 2009 20:06:11 -0000	1.3
+++ mortsrc/common/crdalgo.cpp	27 Sep 2009 16:04:16 -0000
@@ -13,7 +13,7 @@
 
     numvec center( const vector<double>& crd )
     {
-        numvec r(0.0, 3);
+        numvec r(zero_numvec(3));
 
         int natom = crd.size()/3;
         for( int i=0; i < natom; ++i )
@@ -29,7 +29,7 @@
 
     numvec center( const vector<double>& crd, const vector<int>& ids )
     {
-        numvec r(0.0, 6);
+        numvec r(zero_numvec(6));
 
         for( int i=0; i < (int)ids.size(); ++i )
         {
@@ -50,7 +50,8 @@
     {
         assert( crd.size() > 0 );
 
-        numvec r(0.0, 6);
+        numvec r(zero_numvec(6));
+
         r[0] = r[3] = crd[0];
         r[1] = r[4] = crd[1];
         r[2] = r[5] = crd[2];
@@ -73,7 +74,8 @@
     {
         assert( crd.size() > 0 );
 
-        numvec r(0.0, 6);
+        numvec r(zero_numvec(6));
+
         r[0] = DOUBLEMAX; r[3] = -DOUBLEMAX;
         r[1] = DOUBLEMAX; r[4] = -DOUBLEMAX;
         r[2] = DOUBLEMAX; r[5] = -DOUBLEMAX;
@@ -95,7 +97,8 @@
     numvec region( const vector<double>& crd )
     {
         numvec mmx = minmax( crd );
-        numvec rgn(0.0, 6);
+        numvec rgn(zero_numvec(6));
+
         rgn[0] = (mmx[0] + mmx[3])*0.5;
         rgn[1] = (mmx[1] + mmx[4])*0.5;
         rgn[2] = (mmx[2] + mmx[5])*0.5;
@@ -110,7 +113,7 @@
     numvec region( const vector<double>& crd, const vector<double>& vdwr )
     {
         numvec mmx = minmax( crd, vdwr );
-        numvec rgn(0.0,6);
+        numvec rgn(zero_numvec(6));
         rgn[0] = (mmx[0] + mmx[3])*0.5;
         rgn[1] = (mmx[1] + mmx[4])*0.5;
         rgn[2] = (mmx[2] + mmx[5])*0.5;
@@ -262,17 +265,17 @@
     {
         funstack_t::push( "regionlize" );
         numvec rgn = region( crd );
-        translate( crd, -numvec(rgn[slice(0,3,1)]) );
+        translate( crd, -subvec(rgn, 0, 3));
         funstack_t::pop();
-        return rgn[slice(3,3,1)];
+        return subvec(rgn, 3, 3);
     }
 
     numvec regionlize( vector<double>& crd, const vector<double>& vdwr )
     {
         numvec rgn = region( crd, vdwr );
         assert( rgn.size()==6u );
-        translate( crd, -numvec(rgn[slice(0,3,1)]) );
-        return rgn[slice(3,3,1)];
+        translate( crd, -subvec(rgn, 0, 3));
+        return subvec(rgn, 3, 3);
     }
 
     void rotatewald( vector<double>& crd )
Index: mortsrc/common/fortran.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/common/fortran.cpp,v
retrieving revision 1.7
diff -u -r1.7 fortran.cpp
--- mortsrc/common/fortran.cpp	25 Aug 2009 12:52:30 -0000	1.7
+++ mortsrc/common/fortran.cpp	27 Sep 2009 16:04:16 -0000
@@ -129,9 +129,8 @@
 
     void format::write( ostream& os ) const
     {
-        unsigned int bgn = 0;
-        unsigned int vid = 0;
-        int end = m_strfmt.find( '%' );
+        size_t bgn(0), vid(0);
+        size_t end = m_strfmt.find( '%' );
 
         while( end != string::npos )
         {
Index: mortsrc/common/numvec.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/common/numvec.cpp,v
retrieving revision 1.25
diff -u -r1.25 numvec.cpp
--- mortsrc/common/numvec.cpp	16 Feb 2009 03:28:49 -0000	1.25
+++ mortsrc/common/numvec.cpp	27 Sep 2009 16:04:16 -0000
@@ -1,121 +1,103 @@
-#include <iomanip>
-#include <iostream>
+#include <cmath>
 #include <cassert>
-#include <stdexcept>
 #include "fortran.hpp"
 #include "numvec.hpp"
 
-#ifndef M_PI
-#include <math.h>
-#define M_PI 3.1415926535897932385
-#endif
-
-
 namespace mort
 {
-    using std::string;
-
-    using std::runtime_error;
+    numvec makevec(size_t n, const double* x)
+    {
+        numvec result(n);
+        std::copy(x, x + n, result.begin());
+        return result;
+    }
 
-    numvec makevec( double x, double y ) 
+    numvec subvec(const numvec& src, size_t b, size_t n)
     {
-        double r[] = {x, y};
-        return numvec( r, 2 );
+        assert(n > 0);
+        assert(b + n <= src.size());
+        numvec result(n);
+        std::copy(src.begin() + b, src.begin() + b + n, result.begin());
+        return result;
     }
 
-    numvec makevec( double x, double y, double z ) 
+    numvec makevec(const double& x, const double& y)
     {
-        double r[] = {x, y, z};
-        return numvec( r, 3 );
+        numvec result(2);
+        result[0] = x;
+        result[1] = y;
+        return result;
     }
 
-    numvec makevec( double x, double y, double z, double s ) 
+    numvec makevec(const double& x, const double& y, const double& z)
     {
-        double r[] = { x, y, z, s};
-        return numvec( r, 4 );
+        numvec result(3);
+        result[0] = x;
+        result[1] = y;
+        result[2] = z;
+        return result;
     }
 
-    numvec makevec( double x, double y, double z, double s, double u )
+    numvec makevec(const double& x, const double& y, const double& z,
+                   const double& s)
     {
-        double r[] = { x, y, z, s, u };
-        return numvec( r, 5 );
+        numvec result(4);
+        result[0] = x;
+        result[1] = y;
+        result[2] = z;
+        result[3] = s;
+        return result;
     }
 
-    numvec makevec( double x, double y, double z, double s, double u, double v ) 
+    double max(const numvec& v)
     {
-        double r[] = { x, y, z, s, u, v};
-        return numvec( r, 6 );
+        assert(v.size() > 0);
+
+        double accu(v[0]);
+        for (numvec::const_iterator i = v.begin() + 1; i != v.end(); ++i)
+            accu = std::max(accu, *i);
+
+        return accu;
     }
-    
-    double norm( const numvec& v )
+
+    double norm(const numvec& v)
     {
-        if( v.size()==0 )
-            throw runtime_error( "Error: in norm vector size==0" );
+        assert(v.size() > 0);
 
-        return sqrt( dotpd(v,v) );
+        double accu(0);
+        for (numvec::const_iterator i = v.begin(); i != v.end(); ++i)
+            accu += (*i)*(*i);
+
+        return std::sqrt(accu);
     }
 
-    double dotpd( const numvec& a, const numvec& b )
+    double dotpd(const numvec& a, const numvec& b)
     {
-        if( a.size()!=b.size() )
-            throw runtime_error( "Error: dot product two vectors difference size" );
-        
-        return (a*b).sum();
-    }    
- 
-    numvec cross( const numvec& v1, const numvec& v2 )    
+        assert(a.size() == b.size());
+
+        double accu(0);
+        for (size_t n(0); n < a.size(); ++n)
+            accu += a[n]*b[n];
+
+        return accu;
+    }
+
+    numvec cross(const numvec& v1, const numvec& v2)
     {
-        return makevec( v1[1] * v2[2] - v2[1] * v1[2],
-                        v1[2] * v2[0] - v2[2] * v1[0],
-                        v1[0] * v2[1] - v2[0] * v1[1] );
+        return makevec(v1[1]*v2[2] - v2[1]*v1[2],
+                       v1[2]*v2[0] - v2[2]*v1[0],
+                       v1[0]*v2[1] - v2[0]*v1[1]);
     }
 
-    numvec& normalize( numvec& vec )
+    numvec& normalize(numvec& vec)
     {
-        vec /= norm(vec);
-        return vec;
+        return (vec /= norm(vec));
     }
 
-    numvec normalcpy( const numvec& vec )
+    numvec normalcpy(const numvec& vec)
     {
-        numvec copy( vec );
-        normalize( copy );
-        return copy;
+        numvec copy(vec);
+        return normalize(copy);
     }
 
-    
-   
 } // namespace mort
-
-
-
-/*
-  numvec atov( const string& str )
-  {
-  int curly = str.find( '{' );
-
-  int start = (curly == (int)string::npos) ? 0 : curly + 1;
-        
-  std::istringstream is( str.c_str() + start );
-        
-  numvec v;
-  double t;
-  while( is >> t )
-  {
-  v.push_back(t);
-  }
-	
-  return v;
-  }
-
-
-  numvec atov3d( const string& str )
-  {
-  return atov< numvec >( str );
-  }
-    
-  vector4d atov4d( const string& str )
-  {
-  return atov< vector4d >( str );
-  }
-*/
Index: mortsrc/common/numvec.hpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/common/numvec.hpp,v
retrieving revision 1.12
diff -u -r1.12 numvec.hpp
--- mortsrc/common/numvec.hpp	24 Sep 2009 16:59:55 -0000	1.12
+++ mortsrc/common/numvec.hpp	27 Sep 2009 16:04:16 -0000
@@ -1,78 +1,44 @@
 #ifndef MORT_OBJECT_NUMVEC_HPP
 #define MORT_OBJECT_NUMVEC_HPP
 
-#include <vector>
-#include <valarray>
-
+#include <boost/numeric/ublas/vector.hpp>
 
 namespace mort
 {
-    class numvec : public std::valarray<double>
-    {
-
-    public:
-
-        numvec() {}
-
-        numvec( int size ) : std::valarray<double>(size) {}
+    typedef boost::numeric::ublas::vector<double> numvec;
+    typedef boost::numeric::ublas::zero_vector<double> zero_numvec;
+    typedef boost::numeric::ublas::scalar_vector<double> scalar_numvec;
 
-        numvec( double n, int size ) : std::valarray<double>(n, size) {}
+    numvec makevec(size_t, const double*);
 
-        numvec( const double* n, int size ) : std::valarray<double>(n, size) {}
-
-        numvec( const std::valarray<double>& rhs ) : std::valarray<double>( rhs ) {}
-
-        numvec( const std::slice_array<double>& rhs ) :  std::valarray<double>( rhs ) {}
+    /// construct numvec from elements
+    numvec makevec(const double& x, const double& y);
 
-        template< typename T > numvec( const std::_Expr<T,double>& rhs ) : std::valarray<double>( rhs ) {}
+    numvec makevec(const double& x, const double& y, const double& z);
 
-        numvec& operator=( const numvec& rhs ) 
-        { 
-            if( this!=&rhs )
-            {
-                resize( rhs.size() );
-                for( unsigned int i=0; i < size(); ++i )
-                {
-                    (*this)[i] = rhs[i];
-                }
-            }
+    numvec makevec(const double& x, const double& y, const double& z,
+                   const double& s);
 
-            return *this;
-        }
+    numvec subvec(const numvec&, size_t b, size_t n);
 
-        numvec subvec(int bgn, int end) { return (*this)[std::slice(bgn,end,1)]; }
-    };
-  
-    /// construct numvec from elements
-    numvec makevec(double x, double y);
-    
-    numvec makevec(double x, double y, double z);
-
-    numvec makevec(double x, double y, double z, double s);
- 
-    numvec makevec(double x, double y, double z, double s, double u);
-        
-    numvec makevec(double x, double y, double z, double s, double u, double v);        
- 
     // maximum element of a numvec
-    inline double max( const numvec& v ) { return v.max(); }
- 
+    double max(const numvec& vec);
+
     /// norm of a vector 
-    double norm( const numvec& vec );
-  
+    double norm(const numvec& vec);
+
     /// dot product of two vectors.
-    double dotpd( const numvec& v1, const numvec& v2 );
+    double dotpd(const numvec& v1, const numvec& v2);
 
     /// cross product of two vector
-    numvec cross( const numvec& v1, const numvec& v2 );
- 
+    numvec cross(const numvec& v1, const numvec& v2);
+
     /// normalize a vector in place
-    numvec& normalize( numvec& vec );
- 
+    numvec& normalize(numvec& vec);
+
     /// return a normalized copy of a vector, the orginal one remain unchanged.
-    numvec  normalcpy( const numvec& v );    
+    numvec  normalcpy(const numvec& v);
 
 } // namespace mort
 
-#endif
-
+#endif // MORT_OBJECT_NUMVEC_HPP
Index: mortsrc/enefrc/nonbond-pbc.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/enefrc/nonbond-pbc.cpp,v
retrieving revision 1.5
diff -u -r1.5 nonbond-pbc.cpp
--- mortsrc/enefrc/nonbond-pbc.cpp	15 Feb 2009 15:53:39 -0000	1.5
+++ mortsrc/enefrc/nonbond-pbc.cpp	27 Sep 2009 16:04:16 -0000
@@ -119,7 +119,7 @@
 	    }
 	}
 
-	return numvec(edir, eadj);
+	return makevec(edir, eadj);
     }
 
     namespace regewald
Index: mortsrc/object/data.cpp
===================================================================
RCS file: /cvsroot/amber11/src/gleap/mortsrc/object/data.cpp,v
retrieving revision 1.1
diff -u -r1.1 data.cpp
--- mortsrc/object/data.cpp	7 Apr 2009 13:35:27 -0000	1.1
+++ mortsrc/object/data.cpp	27 Sep 2009 16:04:16 -0000
@@ -348,7 +348,7 @@
         {
             map<hashid_t,int>::const_iterator i = m_vlength.find(parmid);
             assert( i != m_vlength.end() );
-            v = numvec(pv, i->second);
+            v = makevec(i->second, pv);
         }
         return result;
     }
