diff --git a/modules/base/doc/generic.rst b/modules/base/doc/generic.rst
index e24ca06b99e5527c73203b3df529c7a81f2a3c61..919497f390775a68deee2fcc973a22f29beace75 100644
--- a/modules/base/doc/generic.rst
+++ b/modules/base/doc/generic.rst
@@ -45,7 +45,7 @@ To store a float value with the key 'myfloatprop' in all atoms of an entity:
   
   import math
   for atom in entity.GetAtomList(): 
-    val=5*math.sin(0.4*atom.GetPos().GetX())
+    val=5*math.sin(0.4*atom.GetPos().x)
     atom.SetFloatProp("myfloatprop", val)
   
 If a GenericProp at a given level (i.e. atom, bond, residue, chain or entity) 
diff --git a/modules/geom/src/composite2.cc b/modules/geom/src/composite2.cc
index b8300fcfb559afe9244583c7b679b1419ef85ead..f6d7f685743c3962d8b8e5acf1a1280b1688b210 100644
--- a/modules/geom/src/composite2.cc
+++ b/modules/geom/src/composite2.cc
@@ -88,16 +88,16 @@ Vec2 Ellipse2::calc_(Real t) const
 Rectangle2 Ellipse2::GetBoundingBox() const
 {
   if(gamma_==0.0 || gamma_==M_PI){
-    return Rectangle2(Vec2(origin_.GetX()-a_,origin_.GetY()-b_),Vec2(origin_.GetX()+a_/2,origin_.GetY()+b_));
+    return Rectangle2(Vec2(origin_.x-a_,origin_.y-b_),Vec2(origin_.x+a_/2,origin_.y+b_));
   }
   if(gamma_==M_PI*0.5 || gamma_==3*M_PI*0.5){
-    return Rectangle2(Vec2(origin_.GetX()-b_,origin_.GetY()-a_),Vec2(origin_.GetX()+b_/2,origin_.GetY()+a_));
+    return Rectangle2(Vec2(origin_.x-b_,origin_.y-a_),Vec2(origin_.x+b_/2,origin_.y+a_));
   }
   Real t=atan(-b_/a_*tan(gamma_));
   Vec2 vx=calc_(t);
   t=atan(b_/a_/tan(gamma_));
   Vec2 vy=calc_(t);
-  return Rectangle2(Vec2(origin_.GetX()-vx.GetX(),origin_.GetY()-vy.GetY()),Vec2(origin_.GetX()+vx.GetX(),origin_.GetY()+vy.GetY()));
+  return Rectangle2(Vec2(origin_.x-vx.x,origin_.y-vy.y),Vec2(origin_.x+vx.x,origin_.y+vy.y));
 }
 
 Rectangle2::Rectangle2():
@@ -112,11 +112,11 @@ bottomright_(bottomright)
 }
 Real Rectangle2::GetWidth() const
 {
-    return bottomright_.GetX()-topleft_.GetX();
+    return bottomright_.x-topleft_.x;
 }
 Real Rectangle2::GetHeight() const
 {
-    return bottomright_.GetY()-topleft_.GetY();
+    return bottomright_.y-topleft_.y;
 }
 Real Rectangle2::GetArea() const
 {
@@ -207,7 +207,7 @@ Real Polygon2::GetArea() const
     unsigned int ip1=(i+1)%nodecount;
     Vec2 n1=GetNode(i);
     Vec2 n2=GetNode(ip1);
-    area+=n1.GetX()*n2.GetY()-n2.GetX()*n1.GetY();
+    area+=n1.x*n2.y-n2.x*n1.y;
   }
   return area/2.0;
 }
@@ -220,8 +220,8 @@ Vec2 Polygon2::GetCentroid() const
     unsigned int ip1=(i+1)%nodecount;
     Vec2 n1=GetNode(i);
     Vec2 n2=GetNode(ip1);
-    cx+=(n1.GetX()+n2.GetX())*(n1.GetX()*n2.GetY()-n2.GetX()*n1.GetY());
-    cy+=(n1.GetY()+n2.GetY())*(n1.GetX()*n2.GetY()-n2.GetX()*n1.GetY());
+    cx+=(n1.x+n2.x)*(n1.x*n2.y-n2.x*n1.y);
+    cy+=(n1.y+n2.y)*(n1.x*n2.y-n2.x*n1.y);
   }
   Real area=GetArea();
   return Vec2(cx/(6.0*area),cy/(6.0*area));
@@ -255,17 +255,17 @@ Rectangle2 Polygon2::GetBoundingBox() const
   {
     return Rectangle2(Vec2(0.0,0.0),Vec2(0.0,0.0));
   }
-  Real x_min=operator[](0).GetX();
+  Real x_min=operator[](0).x;
   Real x_max=x_min;
-  Real y_min=operator[](0).GetY();
+  Real y_min=operator[](0).y;
   Real y_max=y_min;
   unsigned int nodecount=GetNodeCount();
   for(unsigned int i=0;i<nodecount;++i){
     Vec2 n=operator[](i);
-    x_min=std::min<Real>(n.GetX(),x_min);
-    x_max=std::max<Real>(n.GetX(),x_max);
-    y_min=std::min<Real>(n.GetY(),y_min);
-    y_max=std::max<Real>(n.GetY(),y_max);
+    x_min=std::min<Real>(n.x,x_min);
+    x_max=std::max<Real>(n.x,x_max);
+    y_min=std::min<Real>(n.y,y_min);
+    y_max=std::max<Real>(n.y,y_max);
   }
   return Rectangle2(Vec2(x_min,y_min),Vec2(x_max,y_max));
 }
diff --git a/modules/geom/src/composite2_op.cc b/modules/geom/src/composite2_op.cc
index cc660f821eb782bcb6e56a0c4f4dc01002747e6d..5044771a8f5ca2223520692f8a2206c470c05421 100644
--- a/modules/geom/src/composite2_op.cc
+++ b/modules/geom/src/composite2_op.cc
@@ -88,11 +88,11 @@ bool IsInPolygon(const Polygon2& p, const Vec2& v)
   }
   bool c = false;
   Vec2 old=p.GetNode(p.GetNodeCount()-1);
-  bool oldcomp= v.GetY()<old.GetY();
+  bool oldcomp= v.y<old.y;
   for (Polygon2::const_iterator it=p.begin();it!=p.end(); ++it) {
-    bool comp= v.GetY() < (*it).GetY();
+    bool comp= v.y < (*it).y;
     if( comp ^ oldcomp ) {
-      if( (v.GetX() < (old.GetX() - (*it).GetX()) * (v.GetY() - (*it).GetY()) / (old.GetY() - (*it).GetY()) + (*it).GetX()) ){
+      if( (v.x < (old.x - (*it).x) * (v.y - (*it).y) / (old.y - (*it).y) + (*it).x) ){
         c = !c;
       }
     }
@@ -123,7 +123,7 @@ DLLEXPORT bool IsInRectangle(const Rectangle2& r, const Vec2& v)
 {
   Vec2 start=r.GetStart();
   Vec2 end=r.GetEnd();
-  return v.GetX()>=start.GetX() && v.GetX()<=end.GetX() &&  v.GetY()>=start.GetY() && v.GetY()<=end.GetY();
+  return v.x>=start.x && v.x<=end.x &&  v.y>=start.y && v.y<=end.y;
 }
 DLLEXPORT bool IsInCircle(const Circle2& c, const Vec2& v)
 {
diff --git a/modules/geom/src/mat2.cc b/modules/geom/src/mat2.cc
index ac34069dfe055e2dcce97e0e4a5294209a93b5bd..06e196d8cb987f1e4162345a37e40cdcaa764a9a 100644
--- a/modules/geom/src/mat2.cc
+++ b/modules/geom/src/mat2.cc
@@ -65,20 +65,6 @@ bool Mat2::operator==(const Mat2& rhs) const
     data_[1][1] == rhs.data_[1][1];
 }
 
-
-Real& Mat2::operator()(std::size_t r, std::size_t c)
-{
-  if(r>1 || c>1) throw OutOfRangeException();
-  return data_[r][c];
-}
-
-const Real& Mat2::operator()(std::size_t r, std::size_t c) const
-{
-  if(r>1 || c>1) throw OutOfRangeException();
-  return data_[r][c];
-}
-
-
 Mat2& Mat2::operator+=(const Mat2& rhs)
 {
   data_[0][0]+=rhs(0,0);
diff --git a/modules/geom/src/mat2.hh b/modules/geom/src/mat2.hh
index 9231cb7100e3f35541b59608e11c894751448775..bc53416f2caac2772706d50cacb09a76a39effd3 100644
--- a/modules/geom/src/mat2.hh
+++ b/modules/geom/src/mat2.hh
@@ -21,6 +21,7 @@
 
 #include <cstddef> // for size_t
 #include <ostream>
+#include <cassert>
 
 #include <boost/operators.hpp>
 
@@ -63,9 +64,17 @@ public:
   bool operator==(const Mat2& rhs) const;
 
   //! element access
-  Real& operator()(std::size_t r, std::size_t c);
+  Real& operator()(std::size_t r, std::size_t c)
+  {
+    assert(r<=1 && c<=1);
+    return data_[r][c];
+  }
   //! const element access
-  const Real& operator()(std::size_t r, std::size_t c) const;
+  const Real& operator()(std::size_t r, std::size_t c) const
+  {
+    assert(r<=1 && c<=1);
+    return data_[r][c];
+  }
 
   //! addable op
   Mat2& operator+=(const Mat2& rhs);
diff --git a/modules/geom/src/mat3.cc b/modules/geom/src/mat3.cc
index f9492d3d4b4ec7d209de1cde92e31f2d727ceeb9..a202f20a3795fca2a9e36148e9d7f0dd34789900 100644
--- a/modules/geom/src/mat3.cc
+++ b/modules/geom/src/mat3.cc
@@ -24,132 +24,6 @@
 
 namespace geom {
 
-Mat3::Mat3()
-{
-  this->set(1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0);
-}
-
-Mat3::Mat3(Real i00, Real i01, Real i02,Real i10, Real i11, Real i12,Real i20, Real i21, Real i22)
-{
-  this->set(i00,i01,i02,i10,i11,i12,i20,i21,i22);
-}
-
-Mat3::Mat3(const Mat3& m)
-{
-  this->set(m(0,0),m(0,1),m(0,2),m(1,0),m(1,1),m(1,2),m(2,0),m(2,1),m(2,2));
-}
-
-Mat3::Mat3(const Mat2& m)
-{
-  this->set(m(0,0),m(0,1),0.0, m(1,0),m(1,1),0.0, 0.0,0.0,1.0);
-}
-
-Mat3::Mat3(const Real arr[9])
-{
-  this->set(arr[0],arr[1],arr[2],arr[3],arr[4],arr[5],arr[6],arr[7],arr[8]);
-}
-
-Mat3& Mat3::operator=(const Mat3& m)
-{
-  if(&m!=this) {
-    this->set(m(0,0),m(0,1),m(0,2),m(1,0),m(1,1),m(1,2),m(2,0),m(2,1),m(2,2));
-  }
-  return *this;
-}
-
-Mat3 Mat3::Identity() 
-{
-  static Mat3 i(1.0,0.0,0.0,
-                0.0,1.0,0.0,
-                0.0,0.0,1.0);
-  return i;
-}
-
-bool Mat3::operator==(const Mat3& rhs) const
-{
-  return data_[0][0] ==  rhs.data_[0][0] &&
-    data_[1][0] == rhs.data_[1][0] &&
-    data_[2][0] == rhs.data_[2][0] &&
-    data_[0][1] == rhs.data_[0][1] &&
-    data_[1][1] ==  rhs.data_[1][1] &&
-    data_[2][1] ==  rhs.data_[2][1] &&
-    data_[0][2] ==  rhs.data_[0][2] &&
-    data_[1][2] ==  rhs.data_[1][2] &&
-    data_[2][2] ==  rhs.data_[2][2];
-}
-
-Real& Mat3::operator()(std::size_t r, std::size_t c)
-{
-  if(r>2 || c>2) throw OutOfRangeException();
-  return data_[r][c];
-}
-
-const Real& Mat3::operator()(std::size_t r, std::size_t c) const
-{
-  if(r>2 || c>2) throw OutOfRangeException();
-  return data_[r][c];
-}
-
-
-Mat3& Mat3::operator+=(const Mat3& rhs)
-{
-  data_[0][0]+=rhs(0,0);
-  data_[0][1]+=rhs(0,1);
-  data_[0][2]+=rhs(0,2);
-  data_[1][0]+=rhs(1,0);
-  data_[1][1]+=rhs(1,1);
-  data_[1][2]+=rhs(1,2);
-  data_[2][0]+=rhs(2,0);
-  data_[2][1]+=rhs(2,1);
-  data_[2][2]+=rhs(2,2);
-  return *this;
-}
-Mat3& Mat3::operator-=(const Mat3& rhs)
-{
-  data_[0][0]-=rhs(0,0);
-  data_[0][1]-=rhs(0,1);
-  data_[0][2]-=rhs(0,2);
-  data_[1][0]-=rhs(1,0);
-  data_[1][1]-=rhs(1,1);
-  data_[1][2]-=rhs(1,2);
-  data_[2][0]-=rhs(2,0);
-  data_[2][1]-=rhs(2,1);
-  data_[2][2]-=rhs(2,2);
-  return *this;
-}
-Mat3& Mat3::operator*=(const Real d)
-{
-  data_[0][0]*=d;
-  data_[0][1]*=d;
-  data_[0][2]*=d;
-  data_[1][0]*=d;
-  data_[1][1]*=d;
-  data_[1][2]*=d;
-  data_[2][0]*=d;
-  data_[2][1]*=d;
-  data_[2][2]*=d;
-  return *this;
-}
-Mat3& Mat3::operator/=(const Real d)
-{
-  data_[0][0]/=d;
-  data_[0][1]/=d;
-  data_[0][2]/=d;
-  data_[1][0]/=d;
-  data_[1][1]/=d;
-  data_[1][2]/=d;
-  data_[2][0]/=d;
-  data_[2][1]/=d;
-  data_[2][2]/=d;
-  return *this;
-}
-
-Mat3& Mat3::operator*=(const Mat3& m)
-{
-  (*this)=operator*(*this,m);
-  return *this;
-}
-
 geom::Vec3 Mat3::GetCol(int index) const
 {
   return geom::Vec3(data_[0][index], data_[1][index], data_[2][index]);  
@@ -163,12 +37,6 @@ geom::Vec3 Mat3::GetRow(int index) const
 
 ////////////////////////////////////////////
 // private member functions
-void Mat3::set(Real i00, Real i01, Real i02,Real i10, Real i11, Real i12,Real i20, Real i21, Real i22)
-{
-  data_[0][0]=i00; data_[0][1]=i01; data_[0][2]=i02;
-  data_[1][0]=i10; data_[1][1]=i11; data_[1][2]=i12;
-  data_[2][0]=i20; data_[2][1]=i21; data_[2][2]=i22;
-}
 
 std::ostream& operator<<(std::ostream& os, const Mat3& m)
 {
diff --git a/modules/geom/src/mat3.hh b/modules/geom/src/mat3.hh
index f7d90c3f313dff15a1dffb2bd72b5c80fdc4e7ec..15815863328fe4ca1039ced2b5d48250eef952cd 100644
--- a/modules/geom/src/mat3.hh
+++ b/modules/geom/src/mat3.hh
@@ -21,10 +21,12 @@
 
 #include <cstddef> // for size_t
 #include <ostream>
+#include <cassert>
 
 #include <boost/operators.hpp>
 
 #include <ost/geom/module_config.hh>
+#include <ost/geom/mat2.hh>
 
 namespace geom {
 
@@ -36,12 +38,12 @@ class DLLEXPORT_OST_GEOM Mat3:
     private boost::additive1<Mat3>,
     private boost::multiplicative2<Mat3, Real>
 {
-public:
-  static Mat3 Identity();
-  
-  //! Default initialization, identity matrix
-  Mat3();
-
+public:  
+  //! Default initialization, identity matrix  
+  Mat3()
+  {
+    this->set(1.0,0.0,0.0, 0.0,1.0,0.0, 0.0,0.0,1.0);
+  }
   //! In with 9 values in row-major order
   /*!
     row-major order means that the matrix
@@ -52,37 +54,138 @@ public:
 
     is initialized with (a,b,c, d,e,f, g,h,i)
   */
-  Mat3(Real i00, Real i01, Real i02,Real i10, Real i11, Real i12,Real i20, Real i21, Real i22);
-
-  //! Copy ctor
-  Mat3(const Mat3& m);
-
-  // initialization with 2d matrix
-  explicit Mat3(const Mat2& m);
-
-  //! initialization from array
-  explicit Mat3(const Real[9]);
-
-  //! assignement op
-  Mat3& operator=(const Mat3& m);
-
-  //! comparable
-  bool operator==(const Mat3& rhs) const;
+  Mat3(Real i00, Real i01, Real i02,
+       Real i10, Real i11, Real i12,
+       Real i20, Real i21, Real i22)
+  {
+    this->set(i00,i01,i02,i10,i11,i12,i20,i21,i22);
+  }
+
+  Mat3(const Mat3& m)
+  {
+    this->set(m(0,0),m(0,1),m(0,2),m(1,0),m(1,1),m(1,2),m(2,0),m(2,1),m(2,2));
+  }
+
+  Mat3(const Mat2& m)
+  {
+    this->set(m(0, 0), m(0, 1), Real(0.0),
+              m(1, 0), m(1, 1), Real(0.0),
+              Real(0.0), Real(0.0), Real(1.0));
+  }
+  
+  explicit Mat3(const Real arr[9])
+  {
+    this->set(arr[0],arr[1],arr[2],arr[3],arr[4],arr[5],arr[6],arr[7],arr[8]);
+  }  
 
   //! element access
-  Real& operator()(std::size_t r, std::size_t c);
+  Real& operator()(std::size_t r, std::size_t c)
+  {
+    assert(r<=2 && c<=2);
+    return data_[r][c];
+  }
   //! const element access
-  const Real& operator()(std::size_t r, std::size_t c) const;
-
-  //! addable op
-  Mat3& operator+=(const Mat3& rhs);
-  //! subtractable op
-  Mat3& operator-=(const Mat3& rhs);
-
-  Mat3& operator*=(const Real d);
-  Mat3& operator/=(const Real d);
-
-  Mat3& operator*=(const Mat3& m);
+  const Real& operator()(std::size_t r, std::size_t c) const
+  {
+    assert(r<=2 && c<=2);
+    return data_[r][c];
+  }
+
+  Mat3& operator=(const Mat3& m)
+  {
+    if(&m!=this) {
+      this->set(m(0,0),m(0,1),m(0,2),m(1,0),m(1,1),m(1,2),m(2,0),m(2,1),m(2,2));
+    }
+    return *this;
+  }
+
+  static Mat3 Identity() 
+  {
+    static Mat3 i(1.0,0.0,0.0,
+                  0.0,1.0,0.0,
+                  0.0,0.0,1.0);
+    return i;
+  }
+
+  bool operator==(const Mat3& rhs) const
+  {
+    return data_[0][0] ==  rhs.data_[0][0] &&
+      data_[1][0] == rhs.data_[1][0] &&
+      data_[2][0] == rhs.data_[2][0] &&
+      data_[0][1] == rhs.data_[0][1] &&
+      data_[1][1] ==  rhs.data_[1][1] &&
+      data_[2][1] ==  rhs.data_[2][1] &&
+      data_[0][2] ==  rhs.data_[0][2] &&
+      data_[1][2] ==  rhs.data_[1][2] &&
+      data_[2][2] ==  rhs.data_[2][2];
+  }
+
+  Mat3& operator+=(const Mat3& rhs)
+  {
+    data_[0][0]+=rhs(0,0);
+    data_[0][1]+=rhs(0,1);
+    data_[0][2]+=rhs(0,2);
+    data_[1][0]+=rhs(1,0);
+    data_[1][1]+=rhs(1,1);
+    data_[1][2]+=rhs(1,2);
+    data_[2][0]+=rhs(2,0);
+    data_[2][1]+=rhs(2,1);
+    data_[2][2]+=rhs(2,2);
+    return *this;
+  }
+  Mat3& operator-=(const Mat3& rhs)
+  {
+    data_[0][0]-=rhs(0,0);
+    data_[0][1]-=rhs(0,1);
+    data_[0][2]-=rhs(0,2);
+    data_[1][0]-=rhs(1,0);
+    data_[1][1]-=rhs(1,1);
+    data_[1][2]-=rhs(1,2);
+    data_[2][0]-=rhs(2,0);
+    data_[2][1]-=rhs(2,1);
+    data_[2][2]-=rhs(2,2);
+    return *this;
+  }
+  Mat3& operator*=(const Real d)
+  {
+    data_[0][0]*=d;
+    data_[0][1]*=d;
+    data_[0][2]*=d;
+    data_[1][0]*=d;
+    data_[1][1]*=d;
+    data_[1][2]*=d;
+    data_[2][0]*=d;
+    data_[2][1]*=d;
+    data_[2][2]*=d;
+    return *this;
+  }
+  Mat3& operator/=(const Real d)
+  {
+    data_[0][0]/=d;
+    data_[0][1]/=d;
+    data_[0][2]/=d;
+    data_[1][0]/=d;
+    data_[1][1]/=d;
+    data_[1][2]/=d;
+    data_[2][0]/=d;
+    data_[2][1]/=d;
+    data_[2][2]/=d;
+    return *this;
+  }
+
+  Mat3& operator*=(const Mat3& m)
+  {
+    (*this)=Mat3((*this)(0,0)*m(0,0)+(*this)(0,1)*m(1,0)+(*this)(0,2)*m(2,0),
+                 (*this)(0,0)*m(0,1)+(*this)(0,1)*m(1,1)+(*this)(0,2)*m(2,1),
+                 (*this)(0,0)*m(0,2)+(*this)(0,1)*m(1,2)+(*this)(0,2)*m(2,2),
+                 (*this)(1,0)*m(0,0)+(*this)(1,1)*m(1,0)+(*this)(1,2)*m(2,0),
+                 (*this)(1,0)*m(0,1)+(*this)(1,1)*m(1,1)+(*this)(1,2)*m(2,1),
+                 (*this)(1,0)*m(0,2)+(*this)(1,1)*m(1,2)+(*this)(1,2)*m(2,2),                                   
+                 (*this)(2,0)*m(0,0)+(*this)(2,1)*m(1,0)+(*this)(2,2)*m(2,0),
+                 (*this)(2,0)*m(0,1)+(*this)(2,1)*m(1,1)+(*this)(2,2)*m(2,1),
+                 (*this)(2,0)*m(0,2)+(*this)(2,1)*m(1,2)+(*this)(2,2)*m(2,2));
+    return *this;
+  }
 
   Real* Data() {return &data_[0][0];}
   const Real* Data() const {return &data_[0][0];}
@@ -92,9 +195,16 @@ public:
 private:
   Real data_[3][3];
 
-  void set(Real i00, Real i01, Real i02,Real i10, Real i11, Real i12,Real i20, Real i21, Real i22);
+  void set(Real i00, Real i01, Real i02,
+           Real i10, Real i11, Real i12,
+           Real i20, Real i21, Real i22)
+  {
+    data_[0][0]=i00; data_[0][1]=i01; data_[0][2]=i02;
+    data_[1][0]=i10; data_[1][1]=i11; data_[1][2]=i12;
+    data_[2][0]=i20; data_[2][1]=i21; data_[2][2]=i22;
+  }
 };
-
+  
 DLLEXPORT_OST_GEOM std::ostream& operator<<(std::ostream& o, const Mat3& m);
 
 } // ns geom
diff --git a/modules/geom/src/quat.cc b/modules/geom/src/quat.cc
index 3d5f3f2f9c57477a32b718d474d99f27c07e8546..128d540f097fef6d47531b4593387354ae6d81a7 100644
--- a/modules/geom/src/quat.cc
+++ b/modules/geom/src/quat.cc
@@ -462,21 +462,21 @@ Vec3 Quat::Rotate(const Vec3& vec) const {
 Quat Grassmann(const Quat& lhs, const Quat& rhs)
 {
   return Quat(lhs.GetAngle()*rhs.GetAngle()-
-              lhs.GetAxis().GetX()*rhs.GetAxis().GetX()-
-              lhs.GetAxis().GetY()*rhs.GetAxis().GetY()-
-              lhs.GetAxis().GetZ()*rhs.GetAxis().GetZ(),
-                    lhs.GetAngle()*rhs.GetAxis().GetX()+
-                    lhs.GetAxis().GetX()*rhs.GetAngle()+
-                    lhs.GetAxis().GetY()*rhs.GetAxis().GetZ()-
-                    lhs.GetAxis().GetZ()*rhs.GetAxis().GetY(),
-               lhs.GetAngle()*rhs.GetAxis().GetY()-
-               lhs.GetAxis().GetX()*rhs.GetAxis().GetZ()+
-               lhs.GetAxis().GetY()*rhs.GetAngle()+
-               lhs.GetAxis().GetZ()*rhs.GetAxis().GetX(),
-                    lhs.GetAngle()*rhs.GetAxis().GetZ()+
-                    lhs.GetAxis().GetX()*rhs.GetAxis().GetY()-
-                    lhs.GetAxis().GetY()*rhs.GetAxis().GetX()+
-                    lhs.GetAxis().GetZ()*rhs.GetAngle());
+              lhs.GetAxis().x*rhs.GetAxis().x-
+              lhs.GetAxis().y*rhs.GetAxis().y-
+              lhs.GetAxis().z*rhs.GetAxis().z,
+                    lhs.GetAngle()*rhs.GetAxis().x+
+                    lhs.GetAxis().x*rhs.GetAngle()+
+                    lhs.GetAxis().y*rhs.GetAxis().z-
+                    lhs.GetAxis().z*rhs.GetAxis().y,
+               lhs.GetAngle()*rhs.GetAxis().y-
+               lhs.GetAxis().x*rhs.GetAxis().z+
+               lhs.GetAxis().y*rhs.GetAngle()+
+               lhs.GetAxis().z*rhs.GetAxis().x,
+                    lhs.GetAngle()*rhs.GetAxis().z+
+                    lhs.GetAxis().x*rhs.GetAxis().y-
+                    lhs.GetAxis().y*rhs.GetAxis().x+
+                    lhs.GetAxis().z*rhs.GetAngle());
 }
 
 std::ostream& operator<<(std::ostream& str, const Quat& q)
diff --git a/modules/geom/src/vec2.cc b/modules/geom/src/vec2.cc
index 129a2745c25a2cfcec1dcd0a497a9016858b1dfb..4c30939b555d184824b2a4ead68e1b3e49f7c0ea 100644
--- a/modules/geom/src/vec2.cc
+++ b/modules/geom/src/vec2.cc
@@ -22,158 +22,3 @@
 
 #include "exc.hh"
 
-namespace geom {
-
-Vec2::Vec2()
-{
-  this->set(0.0,0.0);
-}
-
-Vec2::Vec2(Real x, Real y)
-{
-  this->set(x,y);
-}
-
-Vec2::Vec2(const Vec2& v)
-{
-  this->set(v.data_[0],v.data_[1]);
-}
-
-Vec2::Vec2(const Vec3& v)
-{
-  this->set(v[0],v[1]); // z is silenly swallowed
-}
-
-Vec2::Vec2(const Vec4& v) 
-{
-  if(v[3]==0.0) {
-    throw DivideByZeroException("zero w component in Vec4->Vec2 conversion");
-  }
-  this->set(v[0]/v[3],v[1]/v[3]); // z is silently swallowed
-}
-
-Vec2::Vec2(const Real d[2])
-{
-  this->set(d[0],d[1]);
-}
-
-/*Vec2::Vec2(const float d[2])
-{
-  this->set(d[0],d[1]);
-}
-*/
-Vec2& Vec2::operator=(const Vec2& v)
-{
-  if(&v!=this) {
-    this->set(v[0],v[1]);
-  }
-  return *this;
-}
-
-bool Vec2::operator==(const Vec2& v) const
-{
-  return data_[0]==v.data_[0] && 
-    data_[1]==v.data_[1];
-}
-
-
-Real& Vec2::operator[](std::size_t indx)
-{
-  if(indx>1) throw OutOfRangeException();
-  return data_[indx];
-}
-
-
-const Real& Vec2::operator[](std::size_t indx) const
-{
-  if(indx>1) throw OutOfRangeException();
-  return data_[indx];
-}
-
-
-Vec2 Vec2::operator-() const
-{
-  Vec2 nrvo(-data_[0],-data_[1]);
-  return nrvo;
-}
-
-Vec2& Vec2::operator+=(const Vec2& rhs)
-{
-  data_[0]+=rhs[0]; data_[1]+=rhs[1];
-  return *this;
-}
-
-Vec2& Vec2::operator-=(const Vec2& rhs)
-{
-  data_[0]-=rhs[0]; data_[1]-=rhs[1];
-  return *this;
-}
-
-Vec2& Vec2::operator*=(Real d)
-{
-  data_[0]*=d; data_[1]*=d;
-  return *this;
-}
-
-Vec2& Vec2::operator+=(Real d)
-{
-  data_[0]+=d; data_[1]+=d;
-  return *this;
-}
-
-Vec2& Vec2::operator-=(Real d)
-{
-  data_[0]-=d; data_[1]-=d;
-  return *this;
-}
-
-Vec2& Vec2::operator/=(Real d)
-{
-  data_[0]/=d; data_[1]/=d;
-  return *this;
-}
-
-////////////////////////////////////////////
-// private member functions
-
-void Vec2::set(Real x, Real y)
-{
-  data_[0]=x; data_[1]=y;
-}
-
-////////////////////////////////////////////
-// functions
-
-Vec2 operator/(Real d, const Vec2& v)
-{
-  Vec2 nrvo(d/v[0],d/v[1]);
-  return nrvo;
-}
-
-std::ostream& operator<<(std::ostream& os, const Vec2& v)
-{
-  os << "(" << v[0] << "," << v[1] << ")";
-  return os;
-}
-
-Real Vec2::GetX() const
-{
-  return data_[0];
-}
-
-Real Vec2::GetY() const
-{
-  return data_[1];
-}
-
-void Vec2::SetX(Real v)
-{
-  data_[0]=v;
-}
-void Vec2::SetY(Real v)
-{
-  data_[1]=v;
-}
-
-
-} // ns
diff --git a/modules/geom/src/vec2.hh b/modules/geom/src/vec2.hh
index 199b6a7605db622828e8b5616cfc074d626529a9..391e073ff17909e411a3d7eb676a2d7d6026a774 100644
--- a/modules/geom/src/vec2.hh
+++ b/modules/geom/src/vec2.hh
@@ -24,8 +24,9 @@
 #include <vector>
 #include <boost/operators.hpp>
 
-#include <ost/geom/module_config.hh>
+
 #include <ost/config.hh>
+#include <ost/geom/module_config.hh>
 
 namespace geom {
 
@@ -36,7 +37,7 @@ class Vec4;
 /*
   Two dimensional vector class, using Real precision.
 */
-class DLLEXPORT_OST_GEOM Vec2:
+class DLLEXPORT Vec2:
     private boost::equality_comparable<Vec2>,
     private boost::additive<Vec2>,
     private boost::additive<Vec2, Real>,
@@ -44,13 +45,13 @@ class DLLEXPORT_OST_GEOM Vec2:
 {
 public:
   //! Default initialization, all components are set to zero
-  Vec2();
+  Vec2(): x(0), y(0) { }
 
   //! Initialization with x, y and z component
-  Vec2(Real x, Real y);
+  Vec2(Real px, Real py): x(px), y(py) { }
 
   //! copy ctor
-  Vec2(const Vec2& v);
+  Vec2(const Vec2& v): x(v.x), y(v.y) { }
 
   //! explicit initialization with 3D vector
   explicit Vec2(const Vec3& v);
@@ -59,53 +60,120 @@ public:
   explicit Vec2(const Vec4& v);
 
   //! explicit initialization with an array of doubles
-  explicit Vec2(const Real[2]);
+  explicit Vec2(const float v[2]): x(v[0]), y(v[1]) { }
+
 
-#if OST_DOUBLE_PRECISION
   //! explicit initialization with an array of floats
-  explicit Vec2(const float[2]);
-#endif
-  //! assignement op
-  Vec2& operator=(const Vec2& v);
+  explicit Vec2(const double v[2]): x(v[0]), y(v[1]) { }
 
-  //! element access
-  Real& operator[](std::size_t indx);
-  //! const element access
-  const Real& operator[](std::size_t indx) const;
-  //! element access
-  Real GetX() const;
-  Real GetY() const;
-  void SetX(Real v);
-  void SetY(Real v);
+
+  Real GetX() const { return x; }
+  Real GetY() const { return y; }
+
+  void SetX(Real d) { x=d; }
+  void SetY(Real d) { y=d; }
 
   //! comparable
-  bool operator==(const Vec2& rhs) const;
+  bool operator==(const Vec2& rhs) const
+  {
+    return x==rhs.x && y==rhs.y;
+  }
 
+  //! element access
+  Real& operator[](std::size_t indx)
+  {
+    return (&x)[indx];
+  }
+  
+  //! const element access
+  const Real& operator[](std::size_t indx) const
+  {
+    
+    return (&x)[indx];
+  }
+  
   //! addable op
-  Vec2& operator+=(const Vec2& rhs);
-  Vec2& operator+=(Real d);
+  Vec2& operator+=(const Vec2& rhs)
+  {
+    x+=rhs.x;
+    y+=rhs.y;
+    return *this;
+  }
+  
+  Vec2& operator+=(Real d)
+  {
+    x+=d;
+    y+=d;
+    return *this;
+  }
+  
   //! subtractable op
-  Vec2& operator-=(const Vec2& rhs);
-  Vec2& operator-=(Real d);
+  Vec2& operator-=(const Vec2& rhs)
+  {
+    x-=rhs.x;
+    y-=rhs.y;
+    return *this;
+  }
+  
+  Vec2& operator-=(Real d)
+  {
+    x-=d;
+    y-=d;
+    return *this;
+  }
   //! negateable
-  Vec2 operator-() const;
+  Vec2 operator-() const
+  {
+    return Vec2(-x, -y);
+  }
+
   //! multipliable
-  Vec2& operator*=(Real d);
+  Vec2& operator*=(Real d)
+  {
+    x*=d;
+    y*=d;
+    return *this;
+  }
+
   //! dividable
-  Vec2& operator/=(Real d);
+  Vec2& operator/=(Real d)
+  {
+    Real one_over_d=Real(1.0)/d;
+    x*=one_over_d;
+    y*=one_over_d;
+    return *this;
+  }
+
+  Real* Data() {return &x;}
+  const Real* Data() const { return &x; }
+
+  Real x;
+  Real y;
+};
 
-  Real* Data() {return data_;}
-  const Real* Data() const {return data_;}
+inline Vec2 operator/(Real d, const Vec2& v)
+{
+  return Vec2(d/v.x, d/v.y);
+}
 
-private:
-  Real data_[2];
+inline std::ostream& operator<<(std::ostream& os, const Vec2& v)
+{
+  os << "(" << v[0] << "," << v[1] << ")";
+  return os;  
+}
+
+}
+
+#include <ost/geom/vec3.hh>
+#include <ost/geom/vec4.hh>
+
+namespace geom {
+
+inline Vec2::Vec2(const Vec3& v): x(v.x), y(v.y) { }
 
-  void set(Real x, Real y);
-};
 
-DLLEXPORT_OST_GEOM Vec2 operator/(Real d, const Vec2& v);
+inline Vec2::Vec2(const Vec4& v): x(v.x), y(v.y) { }
 
-DLLEXPORT_OST_GEOM std::ostream& operator<<(std::ostream&, const Vec2&);
 
 typedef std::vector<Vec2> Vec2List;
 
diff --git a/modules/geom/src/vec3.cc b/modules/geom/src/vec3.cc
index e84c2b64a73f5bfe4a4d712028b1454127474f59..8df7803ce70940310a32e7fe7151419d8ec54691 100644
--- a/modules/geom/src/vec3.cc
+++ b/modules/geom/src/vec3.cc
@@ -22,181 +22,4 @@
 
 #include "exc.hh"
 #include <ost/config.hh>
-namespace geom {
-
-Vec3::Vec3()
-{
-  this->set(0.0,0.0,0.0);
-}
-
-Vec3::Vec3(Real x, Real y, Real z)
-{
-  this->set(x,y,z);
-}
-
-Vec3::Vec3(Real v)
-{
-  this->set(v, v, v);
-}
-
-
-//#if OST_DOUBLE_PRECISION
-//Vec3::Vec3(float x, float y, float z)
-//{
-//  this->set(x,y,z);
-//}
-//#endif
-
-Vec3::Vec3(const Vec3& v)
-{
-  this->set(v.data_[0],v.data_[1],v.data_[2]);
-}
-
-Vec3::Vec3(const Vec2& v)
-{
-  this->set(v[0],v[1],0.0);
-}
-
-Vec3::Vec3(const Vec4& v)
-{
-  if(v[3]==0.0) {
-    throw DivideByZeroException("zero w component in Vec4->Vec3 conversion");
-  }
-  this->set(v[0]/v[3],v[1]/v[3],v[2]/v[3]);
-}
-
-
-Vec3::Vec3(const Real d[3])
-{
-  this->set(d[0],d[1],d[2]);
-}
-#if OST_DOUBLE_PRECISION
-Vec3::Vec3(const float d[3])
-{
-  this->set(d[0],d[1],d[2]);
-}
-#endif
-
-Vec3& Vec3::operator=(const Vec3& v)
-{
-  if(&v!=this) {
-    this->set(v[0],v[1],v[2]);
-  }
-  return *this;
-}
-
-bool Vec3::operator==(const Vec3& v) const
-{
-  return data_[0]==v.data_[0] &&
-    data_[1]==v.data_[1] &&
-    data_[2]==v.data_[2];
-}
-
-Real& Vec3::operator[](std::size_t indx)
-{
-  if(indx>2) throw OutOfRangeException();
-  return data_[indx];
-}
-
-
-const Real& Vec3::operator[](std::size_t indx) const
-{
-  if(indx>2) throw OutOfRangeException();
-  return data_[indx];
-}
-
-Real Vec3::GetX() const
-{
-  return data_[0];
-}
-
-Real Vec3::GetY() const
-{
-  return data_[1];
-}
-
-Real Vec3::GetZ() const
-{
-  return data_[2];
-}
-
-void Vec3::SetX(Real v)
-{
-  data_[0]=v;
-}
-void Vec3::SetY(Real v)
-{
-  data_[1]=v;
-}
-void Vec3::SetZ(Real v)
-{
-  data_[2]=v;
-}
-
-Vec3 Vec3::operator-() const
-{
-  Vec3 nrvo(-data_[0],-data_[1],-data_[2]);
-  return nrvo;
-}
-
-Vec3& Vec3::operator+=(const Vec3& rhs)
-{
-  data_[0]+=rhs[0]; data_[1]+=rhs[1]; data_[2]+=rhs[2];
-  return *this;
-}
-
-Vec3& Vec3::operator-=(const Vec3& rhs)
-{
-  data_[0]-=rhs[0]; data_[1]-=rhs[1]; data_[2]-=rhs[2];
-  return *this;
-}
-
-Vec3& Vec3::operator*=(Real d)
-{
-  data_[0]*=d; data_[1]*=d; data_[2]*=d;
-  return *this;
-}
-
-Vec3& Vec3::operator+=(Real d)
-{
-  data_[0]+=d; data_[1]+=d; data_[2]+=d;
-  return *this;
-}
-
-Vec3& Vec3::operator-=(Real d)
-{
-  data_[0]-=d; data_[1]-=d; data_[2]-=d;
-  return *this;
-}
-
-Vec3& Vec3::operator/=(Real d)
-{
-  data_[0]/=d; data_[1]/=d; data_[2]/=d;
-  return *this;
-}
-
-////////////////////////////////////////////
-// private member functions
-
-void Vec3::set(Real x, Real y, Real z)
-{
-  data_[0]=x; data_[1]=y; data_[2]=z;
-}
-
-////////////////////////////////////////////
-// functions
-
-Vec3 operator/(Real d, const Vec3& v)
-{
-  Vec3 nrvo(d/v[0],d/v[1],d/v[2]);
-  return nrvo;
-}
-
-std::ostream& operator<<(std::ostream& os, const Vec3& v)
-{
-  os << "(" << v[0] << "," << v[1] << "," << v[2] << ")";
-  return os;
-}
-
-} // ns geom
 
diff --git a/modules/geom/src/vec3.hh b/modules/geom/src/vec3.hh
index ac64e76f780564110bb1e34b04ab1174c7495d9b..627de4c1ee9cea9f6da151cd98a1a70879ab8c19 100644
--- a/modules/geom/src/vec3.hh
+++ b/modules/geom/src/vec3.hh
@@ -24,8 +24,9 @@
 #include <vector>
 #include <boost/operators.hpp>
 
-#include <ost/geom/module_config.hh>
+
 #include <ost/config.hh>
+#include <ost/geom/module_config.hh>
 
 namespace geom {
 
@@ -43,13 +44,13 @@ class DLLEXPORT_OST_GEOM Vec3:
 {
 public:
   //! Default initialization, all components are set to zero
-  Vec3();
+  Vec3(): x(0), y(0), z(0) {}
 
   //! Initialization with x, y and z component
-  Vec3(Real x, Real y, Real z);
+  Vec3(Real px, Real py, Real pz): x(px), y(py), z(pz) {}
 
   //! copy ctor
-  Vec3(const Vec3& v);
+  Vec3(const Vec3& v): x(v.x), y(v.y), z(v.z) { }
 
   //! (implicit) initialization with 2D vector
   Vec3(const Vec2& v);
@@ -60,62 +61,139 @@ public:
     to a 3D vector, resulting in (x/w,y/w,z/w)
   */
   explicit Vec3(const Vec4& v);
-
+  
+  explicit Vec3(Real v): x(v), y(v), z(v) { }
+  
   //! explicit initialization with an array of doubles
-  explicit Vec3(const Real[3]);
+  explicit Vec3(const double v[3]): x(v[0]), y(v[1]), z(v[2]) { }
 
-  explicit Vec3(Real v);
-  
-#if OST_DOUBLE_PRECISION
   //! explicit initialization with an array of floats
-  explicit Vec3(const float[3]);
-#endif
+  explicit Vec3(const float v[3]): x(v[0]), y(v[1]), z(v[2]) { }
+
   //! assignement op
-  Vec3& operator=(const Vec3& v);
+  Vec3& operator=(const Vec3& v)
+  {
+    x=v.x;
+    y=v.y;
+    z=v.z;
+    return *this;
+  }
 
   //! comparable
-  bool operator==(const Vec3& rhs) const;
+  bool operator==(const Vec3& rhs) const
+  {
+    return x==rhs.x && y==rhs.y && z==rhs.z;
+  }
 
   //! element access
-  Real& operator[](std::size_t indx);
+  Real& operator[](std::size_t indx)
+  {
+    return (&x)[indx];
+  }
+  
   //! const element access
-  const Real& operator[](std::size_t indx) const;
+  const Real& operator[](std::size_t indx) const
+  {
+    
+    return (&x)[indx];
+  }
   //! element access
-  Real GetX() const;
-  Real GetY() const;
-  Real GetZ() const;
-  void SetX(Real v);
-  void SetY(Real v);
-  void SetZ(Real v);
+  Real GetX() const { return x; }
+  Real GetY() const { return y; }
+  Real GetZ() const { return z; }
+  void SetX(Real v) { x=v; }
+  void SetY(Real v) { y=v; }
+  void SetZ(Real v) { z=v; }
 
   //! addable op
-  Vec3& operator+=(const Vec3& rhs);
-  Vec3& operator+=(Real d);
+  Vec3& operator+=(const Vec3& rhs)
+  {
+    x+=rhs.x;
+    y+=rhs.y;
+    z+=rhs.z;
+    return *this;
+  }
+  
+  Vec3& operator+=(Real d)
+  {
+    x+=d;
+    y+=d;
+    z+=d;
+    return *this;
+  }
+  
   //! subtractable op
-  Vec3& operator-=(const Vec3& rhs);
-  Vec3& operator-=(Real d);
+  Vec3& operator-=(const Vec3& rhs)
+  {
+    x-=rhs.x;
+    y-=rhs.y;
+    z-=rhs.z;
+    return *this;
+  }
+  
+  Vec3& operator-=(Real d)
+  {
+    x-=d;
+    y-=d;
+    z-=d;
+    return *this;
+  }
   //! negateable
-  Vec3 operator-() const;
+  Vec3 operator-() const
+  {
+    return Vec3(-x, -y, -z);
+  }
+  
   //! multipliable
-  Vec3& operator*=(Real d);
+  Vec3& operator*=(Real d)
+  {
+    x*=d;
+    y*=d;
+    z*=d;
+    return *this;
+  }
+  
   //! dividable
-  Vec3& operator/=(Real d);
-
-  Real* Data() {return data_;}
-  const Real* Data() const {return data_;}
+  Vec3& operator/=(Real d)
+  {
+    Real one_over_d=Real(1.0)/d;
+    x*=one_over_d;
+    y*=one_over_d;
+    z*=one_over_d;
+    return *this;
+  }
+
+  Real* Data() {return &x;}
+  const Real* Data() const {return &x;}
+
+  Real x;
+  Real y;
+  Real z;
+};
 
-private:
-  Real data_[3];
+inline Vec3 operator/(Real d, const Vec3& v)
+{
+  Vec3 nrvo(d/v[0],d/v[1],d/v[2]);
+  return nrvo;
+}
 
-  void set(Real x, Real y, Real z);
-};
+inline std::ostream& operator<<(std::ostream& os, const Vec3& v)
+{
+  os << "[" << v.x << ", " << v.y << "]";
+  return os;
+}
+}
 
-DLLEXPORT_OST_GEOM Vec3 operator/(Real, const Vec3& v);
+#include <ost/geom/vec2.hh>
+#include <ost/geom/vec4.hh>
 
-DLLEXPORT_OST_GEOM std::ostream& operator<<(std::ostream&, const Vec3&);
+namespace geom {
 
 typedef std::vector<Vec3> Vec3List;
 
+inline Vec3::Vec3(const Vec2& v): x(v.x), y(v.y), z(0.0) { }
+inline Vec3::Vec3(const Vec4& v): x(v.x/v.w), y(v.y/v.w), z(v.z/v.w) { }
+  
 } // namespace geom
 
 
diff --git a/modules/geom/src/vec4.cc b/modules/geom/src/vec4.cc
index 8f53869a97983bea6f0ae6adc116327cfe9323bd..4c30939b555d184824b2a4ead68e1b3e49f7c0ea 100644
--- a/modules/geom/src/vec4.cc
+++ b/modules/geom/src/vec4.cc
@@ -22,136 +22,3 @@
 
 #include "exc.hh"
 
-namespace geom {
-
-Vec4::Vec4()
-{
-  this->set(0.0,0.0,0.0,1.0);
-}
-
-Vec4::Vec4(Real x, Real y, Real z, Real w)
-{
-  this->set(x,y,z,w);
-}
-
-Vec4::Vec4(const Vec4& v)
-{
-  this->set(v.data_[0],v.data_[1],v.data_[2],v.data_[3]);
-}
-
-Vec4::Vec4(const Vec2& v)
-{
-  this->set(v[0],v[1],0.0,1.0);
-}
-
-Vec4::Vec4(const Vec3& v)
-{
-  this->set(v[0],v[1],v[2],1.0);
-}
-
-Vec4::Vec4(const Real d[4])
-{
-  this->set(d[0],d[1],d[2],d[3]);
-}
-
-#if OST_DOUBLE_PRECISION
-Vec4::Vec4(const float d[4])
-{
-  this->set(d[0],d[1],d[2],d[3]);
-}
-#endif
-
-Vec4& Vec4::operator=(const Vec4& v)
-{
-  if(&v!=this) {
-    this->set(v[0],v[1],v[2],v[3]);
-  }
-  return *this;
-}
-
-bool Vec4::operator==(const Vec4& v) const
-{
-  return data_[0]==v.data_[0] && 
-    data_[1]==v.data_[1] &&
-    data_[2]==v.data_[2] &&
-    data_[3]==v.data_[3];
-}
-
-Real& Vec4::operator[](std::size_t indx)
-{
-  if(indx>3) throw OutOfRangeException();
-  return data_[indx];
-}
-
-const Real& Vec4::operator[](std::size_t indx) const
-{
-  if(indx>3) throw OutOfRangeException();
-  return data_[indx];
-}
-
-Vec4 Vec4::operator-() const
-{
-  Vec4 nrvo(-data_[0],-data_[1],-data_[2],-data_[3]);
-  return nrvo;
-}
-
-Vec4& Vec4::operator+=(const Vec4& rhs)
-{
-  data_[0]+=rhs[0]; data_[1]+=rhs[1]; data_[2]+=rhs[2]; data_[3]+=rhs[3];
-  return *this;
-}
-
-Vec4& Vec4::operator-=(const Vec4& rhs)
-{
-  data_[0]-=rhs[0]; data_[1]-=rhs[1]; data_[2]-=rhs[2]; data_[3]-=rhs[3];
-  return *this;
-}
-
-Vec4& Vec4::operator*=(Real d)
-{
-  data_[0]*=d; data_[1]*=d; data_[2]*=d; data_[3]*=d;
-  return *this;
-}
-
-Vec4& Vec4::operator+=(Real d)
-{
-  data_[0]+=d; data_[1]+=d; data_[2]+=d; data_[3]+=d;
-  return *this;
-}
-
-Vec4& Vec4::operator-=(Real d)
-{
-  data_[0]-=d; data_[1]-=d; data_[2]-=d; data_[3]-=d;
-  return *this;
-}
-
-Vec4& Vec4::operator/=(Real d)
-{
-  data_[0]/=d; data_[1]/=d; data_[2]/=d; data_[3]/=d;
-  return *this;
-}
-
-////////////////////////////////////////////
-// private member functions
-
-void Vec4::set(Real x, Real y, Real z, Real w)
-{
-  data_[0]=x; data_[1]=y; data_[2]=z; data_[3]=w;
-}
-
-////////////////////////////////////////////
-// functions
-
-Vec4 operator/(Real d, const Vec4& v)
-{
-  Vec4 nrvo(d/v[0],d/v[1],d/v[2],d/v[3]);
-  return nrvo;
-}
-
-std::ostream& operator<<(std::ostream& os, const Vec4& v)
-{
-  os << "(" << v[0] << "," << v[1] << "," << v[2] << "," << v[3] << ")";
-  return os;
-}
-
-}
diff --git a/modules/geom/src/vec4.hh b/modules/geom/src/vec4.hh
index b114713ddd232a269b892da4edd62044ba202c3a..ad3367b97b8f7f30f95aab6c5da54ccc91cb987d 100644
--- a/modules/geom/src/vec4.hh
+++ b/modules/geom/src/vec4.hh
@@ -32,12 +32,11 @@ namespace geom {
 // fw decl
 class Vec2;
 class Vec3;
-class Mat4;
 
 /*
   Four dimensional, homogeneous vector class, using Real precision.
 */
-class DLLEXPORT_OST_GEOM Vec4:
+class DLLEXPORT Vec4:
     private boost::equality_comparable<Vec4>,
     private boost::additive<Vec4>,
     private boost::additive<Vec4, Real>,
@@ -45,13 +44,13 @@ class DLLEXPORT_OST_GEOM Vec4:
 {
 public:
   //! Default initialization, all components are set to zero
-  Vec4();
+  Vec4(): x(0), y(0), z(0), w(0) { }
 
   //! Initialization with x, y and z component
-  Vec4(Real x, Real y, Real z, Real w);
+  Vec4(Real px, Real py, Real pz, Real pw): x(px), y(py), z(pz), w(pw) { }
 
   //! copy ctor
-  Vec4(const Vec4& v);
+  Vec4(const Vec4& v): x(v[0]), y(v[1]), z(v[2]), w(v[3]) { }
 
   //! (implicit) initialization with 2D vector
   Vec4(const Vec2& v);
@@ -59,50 +58,146 @@ public:
   //! (implicit) initialization with 3D vector
   Vec4(const Vec3& v);
 
-  //! explicit initialization with an array of doubles
-  explicit Vec4(const Real[4]);
-
-#if OST_DOUBLE_PRECISION
   //! explicit initialization with an array of floats
-  explicit Vec4(const float[4]);
-#endif
-
+  explicit Vec4(const float v[4]): x(v[0]), y(v[1]), z(v[2]), w(v[3]) { }
+  
+  //! explicit initialization with an array of doubles  
+  explicit Vec4(const double v[4]): x(v[0]), y(v[1]), z(v[2]), w(v[3]) { }
   //! assignement op
-  Vec4& operator=(const Vec4& v);
+  Vec4& operator=(const Vec4& v)
+  {
+    x=v.x;
+    y=v.y;
+    z=v.z;
+    w=v.w;
+    return *this;
+  }
 
   //! comparable
-  bool operator==(const Vec4& rhs) const;
+  bool operator==(const Vec4& rhs) const
+  {
+    return x==rhs.x && y==rhs.y && z==rhs.z && w==rhs.w;
+  }
 
   //! element access
-  Real& operator[](std::size_t indx);
+  Real& operator[](std::size_t indx)
+  {
+    return (&x)[indx];
+  }
+  
   //! const element access
-  const Real& operator[](std::size_t indx) const;
+  const Real& operator[](std::size_t indx) const
+  {
+    
+    return (&x)[indx];
+  }
+  //! element access
+  Real GetX() const { return x; }
+  Real GetY() const { return y; }
+  Real GetZ() const { return z; }
+  Real GetW() const { return w; }
+  
+  void SetX(Real v) { x=v; }
+  void SetY(Real v) { y=v; }
+  void SetZ(Real v) { z=v; }
+  void SetW(Real v) { w=v; }  
 
   //! addable op
-  Vec4& operator+=(const Vec4& rhs);
-  Vec4& operator+=(Real d);
+  Vec4& operator+=(const Vec4& rhs)
+  {
+    x+=rhs.x;
+    y+=rhs.y;
+    z+=rhs.z;
+    w+=rhs.w;
+    return *this;
+  }
+  
+  Vec4& operator+=(Real d)
+  {
+    x+=d;
+    y+=d;
+    z+=d;
+    w+=d;
+    return *this;
+  }
+  
   //! subtractable op
-  Vec4& operator-=(const Vec4& rhs);
-  Vec4& operator-=(Real d);
+  Vec4& operator-=(const Vec4& rhs)
+  {
+    x-=rhs.x;
+    y-=rhs.y;
+    z-=rhs.z;
+    w-=rhs.w;
+    return *this;
+  }
+  
+  Vec4& operator-=(Real d)
+  {
+    x-=d;
+    y-=d;
+    z-=d;
+    w-=d;
+    return *this;
+  }
   //! negateable
-  Vec4 operator-() const;
+  Vec4 operator-() const
+  {
+    return Vec4(-x, -y, -z, -w);
+  }
+  
   //! multipliable
-  Vec4& operator*=(Real d);
+  Vec4& operator*=(Real d)
+  {
+    x*=d;
+    y*=d;
+    z*=d;
+    w*=d;
+    return *this;
+  }
+  
   //! dividable
-  Vec4& operator/=(Real d);
+  Vec4& operator/=(Real d)
+  {
+    Real one_over_d=Real(1.0)/d;
+    x*=one_over_d;
+    y*=one_over_d;
+    z*=one_over_d;
+    w*=one_over_d;
+    return *this;
+  }
+
+  Real* Data() {return &x;}
+  const Real* Data() const {return &x;}
+
+  Real x;
+  Real y;
+  Real z;
+  Real w;
+};
 
-  Real* Data() {return data_;}
-  const Real* Data() const {return data_;}
+inline Vec4 operator/(Real d, const Vec4& v)
+{
+  Vec4 nrvo(d/v[0],d/v[1],d/v[2],d/v[3]);
+  return nrvo;
+}
 
-private:
-  Real data_[4];
+inline std::ostream& operator<<(std::ostream& os, const Vec4& v)
+{
+  os << "(" << v[0] << "," << v[1] << "," << v[2] << "," << v[3] << ")";
+  return os;
+}
+}
 
-  void set(Real x, Real y, Real z, Real w);
-};
+#include <ost/geom/vec2.hh>
+#include <ost/geom/vec3.hh>
 
-DLLEXPORT_OST_GEOM Vec4 operator/(Real, const Vec4& v);
+namespace geom {
+  
+//! (implicit) initialization with 2D vector
+inline Vec4::Vec4(const Vec2& v): x(v.x), y(v.y), z(0), w(1) { }
 
-DLLEXPORT_OST_GEOM std::ostream& operator<<(std::ostream&, const Vec4&);
+//! (implicit) initialization with 3D vector
+inline Vec4::Vec4(const Vec3& v): x(v.x), y(v.y), z(v.z), w(1.0) { }
 
 } // namespace geom
 
diff --git a/modules/geom/src/vecmat2_op.cc b/modules/geom/src/vecmat2_op.cc
index 6bda1a12a5cdbadc6735c8d2f5e665a91bcc1133..c23866293e72f928ed79d7938ad45c22f0b46d19 100644
--- a/modules/geom/src/vecmat2_op.cc
+++ b/modules/geom/src/vecmat2_op.cc
@@ -26,70 +26,6 @@
 #include "mat2.hh"
 
 namespace geom {
-
-Real Length(const Vec2& v)
-{
-  return std::sqrt(Length2(v));
-}
-
-Real Length2(const Vec2& v)
-{
-  return v[0]*v[0]+v[1]*v[1];
-}
-
-bool Equal(const Vec2& v1, const Vec2& v2, Real ephilon)
-{
-  return std::fabs(v1[0]-v2[0])<ephilon &&
-    std::fabs(v1[1]-v2[1])<ephilon;
-}
-bool Equal(const Mat2& m1, const Mat2& m2, Real ephilon)
-{
-  return std::fabs(m1(0,0)-m2(0,0))<ephilon &&
-    std::fabs(m1(0,1)-m2(0,1))<ephilon &&
-    std::fabs(m1(1,0)-m2(1,0))<ephilon &&
-    std::fabs(m1(1,1)-m2(1,1))<ephilon;
-}
-
-Real Dot(const Vec2& v1, const Vec2& v2)
-{
-  return v1[0]*v2[0]+v1[1]*v2[1];
-}
-
-Vec2 Normalize(const Vec2& v)
-{
-  Real l=Length(v);
-  if(l==0.0) {
-    return v;
-  }
-  return v/l;
-}
-
-
-Vec2 CompMultiply(const Vec2& v1, const Vec2& v2)
-{
-  Vec2 nrvo(v1[0]*v2[0],v1[1]*v2[1]);
-  return nrvo;
-}
-
-Vec2 CompDivide(const Vec2& v1, const Vec2& v2)
-{
-  Vec2 nrvo(v1[0]/v2[0],v1[1]/v2[1]);
-  return nrvo;
-}
-
-Vec2 operator*(const Vec2& v,const Mat2& m)
-{
-  Vec2 nrvo(v[0]*m(0,0)+v[1]*m(1,0),
-      v[0]*m(0,1)+v[1]*m(1,1));
-  return nrvo;
-}
-
-Vec2 operator*(const Mat2& m, const Vec2& v)
-{
-  Vec2 nrvo(v[0]*m(0,0)+v[1]*m(0,1),
-      v[0]*m(1,0)+v[1]*m(1,1));
-  return nrvo;
-}
 // determinant
 Real Det(const Mat2& m)
 {
@@ -116,14 +52,6 @@ Real Angle(const Vec2& v1, const Vec2& v2)
 {
   return std::acos(Dot(v1,v2)/Length(v1)/Length(v2));
 }
-Mat2 operator*(const Mat2& m1, const Mat2& m2)
-{
-  Mat2 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0),
-      m1(0,0)*m2(0,1)+m1(0,1)*m2(1,1),
-      m1(1,0)*m2(0,0)+m1(1,1)*m2(1,0),
-      m1(1,0)*m2(0,1)+m1(1,1)*m2(1,1));
-  return nrvo;
-}
 
 Real SignedAngle(const Vec2& v1, const Vec2& v2)
 {
@@ -131,21 +59,7 @@ Real SignedAngle(const Vec2& v1, const Vec2& v2)
   vc=Cross(Vec3(v1),Vec3(v2));
   if(Length(v1)==0.0 || Length(v2)==0.0 || Length(vc)==0.0)
     return 0.0;
-  return acos(Dot(v1,v2)/Length(v1)/Length(v2))*vc.GetZ()/std::fabs(vc.GetZ());
-}
-
-Vec2 Min(const Vec2& v1, const Vec2& v2)
-{
-  Vec2 nrvo(std::min(v1[0],v2[0]),
-            std::min(v1[1],v2[1]));
-  return nrvo;
-}
-
-Vec2 Max(const Vec2& v1, const Vec2& v2)
-{
-  Vec2 nrvo(std::max(v1[0],v2[0]),
-            std::max(v1[1],v2[1]));
-  return nrvo;
+  return acos(Dot(v1,v2)/Length(v1)/Length(v2))*vc.z/std::fabs(vc.z);
 }
 
 Vec2 Rotate(const Vec2& v,Real ang)
diff --git a/modules/geom/src/vecmat2_op.hh b/modules/geom/src/vecmat2_op.hh
index cf20513678438ae61f390aea724cd8bcc933c6f2..31ae6430aca27c9ade3a3df4ec3a4786fb0c8cdb 100644
--- a/modules/geom/src/vecmat2_op.hh
+++ b/modules/geom/src/vecmat2_op.hh
@@ -23,33 +23,71 @@
 
 #include <ost/geom/module_config.hh>
 
+#include <ost/geom/vec2.hh>
+#include <ost/geom/mat2.hh>
+
 namespace geom {
 
 class Vec2;
 class Mat2;
 
-//! returns length of vector
-Real DLLEXPORT_OST_GEOM Length(const Vec2& v);
-
 //! returns squared length of vector
-Real DLLEXPORT_OST_GEOM Length2(const Vec2& v);
+inline Real Length2(const Vec2& v)
+{
+  return v[0]*v[0]+v[1]*v[1];
+}
+
+//! returns length of vector
+inline Real Length(const Vec2& v)
+{
+  return std::sqrt(Length2(v));
+}
 
 //! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Vec2& v1, const Vec2& v2, Real ephilon=EPSILON);
+inline bool Equal(const Vec2& v1, const Vec2& v2, Real ephilon=EPSILON)
+{
+  return std::fabs(v1[0]-v2[0])<ephilon &&
+    std::fabs(v1[1]-v2[1])<ephilon;
+}
+
 //! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Mat2& m1, const Mat2& m2, Real ephilon=EPSILON);
+inline bool Equal(const Mat2& m1, const Mat2& m2, Real ephilon=EPSILON)
+{
+  return std::fabs(m1(0,0)-m2(0,0))<ephilon &&
+    std::fabs(m1(0,1)-m2(0,1))<ephilon &&
+    std::fabs(m1(1,0)-m2(1,0))<ephilon &&
+    std::fabs(m1(1,1)-m2(1,1))<ephilon;
+}
 
 //! vector dot product
-Real DLLEXPORT_OST_GEOM Dot(const Vec2& v1, const Vec2& v2);
+inline Real Dot(const Vec2& v1, const Vec2& v2)
+{
+    return v1[0]*v2[0]+v1[1]*v2[1];
+}
 
 //! Normalize vector
-Vec2 DLLEXPORT_OST_GEOM Normalize(const Vec2& v);
+inline Vec2 Normalize(const Vec2& v)
+{
+  Real l=Length(v);
+  if(l==0.0) {
+    return v;
+  }
+  return v/l;
+}
 
 //! multiply each component of v1 with that of v2
-Vec2 DLLEXPORT_OST_GEOM CompMultiply(const Vec2& v1, const Vec2& v2);
+inline Vec2 CompMultiply(const Vec2& v1, const Vec2& v2)
+{
+  Vec2 nrvo(v1[0]*v2[0],v1[1]*v2[1]);
+  return nrvo;
+}
 
 //! divide each component of v1 by that of v2
-Vec2 DLLEXPORT_OST_GEOM CompDivide(const Vec2& v1, const Vec2& v2);
+inline Vec2 CompDivide(const Vec2& v1, const Vec2& v2)
+{
+  Vec2 nrvo(v1[0]/v2[0],v1[1]/v2[1]);
+  return nrvo;
+}
 
 //! vector matrix multiplication
 /*!
@@ -58,7 +96,12 @@ Vec2 DLLEXPORT_OST_GEOM CompDivide(const Vec2& v1, const Vec2& v2);
           |a b|
   (x,y) * |c d| = (ax+cy, bx+dy)
 */
-Vec2 DLLEXPORT_OST_GEOM operator*(const Vec2& v,const Mat2& m);
+inline Vec2 operator*(const Vec2& v,const Mat2& m)
+{
+  Vec2 nrvo(v[0]*m(0,0)+v[1]*m(1,0),
+      v[0]*m(0,1)+v[1]*m(1,1));
+  return nrvo;
+}
 
 //! vector matrix multiplication
 /*!
@@ -67,7 +110,12 @@ Vec2 DLLEXPORT_OST_GEOM operator*(const Vec2& v,const Mat2& m);
    |a b|   (x)   (ax+by)
    |c d| * (y)  =(cx+dy)
 */
-Vec2 DLLEXPORT_OST_GEOM operator*(const Mat2& m, const Vec2& v);
+inline Vec2 operator*(const Mat2& m, const Vec2& v)
+{
+  Vec2 nrvo(v[0]*m(0,0)+v[1]*m(0,1),
+      v[0]*m(1,0)+v[1]*m(1,1));
+  return nrvo;
+}
 
 //! determinant
 Real DLLEXPORT_OST_GEOM Det(const Mat2& m);
@@ -83,13 +131,28 @@ Real DLLEXPORT_OST_GEOM Angle(const Vec2& v1, const Vec2& v2);
 Real  DLLEXPORT_OST_GEOM SignedAngle(const Vec2& v1, const Vec2& v2);
 
 //! matrix matrix multiplication
-Mat2 DLLEXPORT_OST_GEOM operator*(const Mat2& m1, const Mat2& m2);
-
-//! returns std::min of each component
-Vec2 DLLEXPORT_OST_GEOM Min(const Vec2& v1, const Vec2& v2);
-
-//! returns std::max of each component
-Vec2 DLLEXPORT_OST_GEOM Max(const Vec2& v1, const Vec2& v2);
+inline Mat2 operator*(const Mat2& m1, const Mat2& m2)
+{
+  Mat2 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0),
+      m1(0,0)*m2(0,1)+m1(0,1)*m2(1,1),
+      m1(1,0)*m2(0,0)+m1(1,1)*m2(1,0),
+      m1(1,0)*m2(0,1)+m1(1,1)*m2(1,1));
+  return nrvo;
+}
+
+inline Vec2 Min(const Vec2& v1, const Vec2& v2)
+{
+  Vec2 nrvo(std::min(v1[0],v2[0]),
+            std::min(v1[1],v2[1]));
+  return nrvo;
+}
+
+inline Vec2 Max(const Vec2& v1, const Vec2& v2)
+{
+  Vec2 nrvo(std::max(v1[0],v2[0]),
+            std::max(v1[1],v2[1]));
+  return nrvo;
+}
 
 //! Rotation
 DLLEXPORT Vec2 Rotate(const Vec2& v,Real ang);
diff --git a/modules/geom/src/vecmat3_op.cc b/modules/geom/src/vecmat3_op.cc
index 5487c5d57bf406eff9e29232fd8583815b5e14b6..8ef740e3769afc4d967696d1a0e88eca11c9ad08 100644
--- a/modules/geom/src/vecmat3_op.cc
+++ b/modules/geom/src/vecmat3_op.cc
@@ -29,102 +29,6 @@
 
 namespace geom {
 
-Real Length(const Vec3& v)
-{
-  return std::sqrt(Length2(v));
-}
-
-Real Length2(const Vec3& v)
-{
-  return v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
-}
-
-bool Equal(const Vec3& v1, const Vec3& v2, Real ephilon)
-{
-  return std::fabs(v1[0]-v2[0])<ephilon &&
-    std::fabs(v1[1]-v2[1])<ephilon &&
-    std::fabs(v1[2]-v2[2])<ephilon;
-}
-
-bool Equal(const Mat3& m1, const Mat3& m2, Real ephilon)
-{
-  return std::fabs(m1(0,0)-m2(0,0))<ephilon &&
-    std::fabs(m1(0,1)-m2(0,1))<ephilon &&
-    std::fabs(m1(0,2)-m2(0,2))<ephilon &&
-    std::fabs(m1(1,0)-m2(1,0))<ephilon &&
-    std::fabs(m1(1,1)-m2(1,1))<ephilon &&
-    std::fabs(m1(1,2)-m2(1,2))<ephilon &&
-    std::fabs(m1(2,0)-m2(2,0))<ephilon &&
-    std::fabs(m1(2,1)-m2(2,1))<ephilon &&
-    std::fabs(m1(2,2)-m2(2,2))<ephilon;
-}
-
-Real Dot(const Vec3& v1, const Vec3& v2)
-{
-  return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
-}
-
-Vec3 Cross(const Vec3& v1, const Vec3& v2)
-{
-  Vec3 nrvo(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 nrvo;
-}
-
-Vec3 Normalize(const Vec3& v)
-{
-  Real l=Length(v);
-  if(l==0.0) {
-    return v;
-  }
-  return v/l;
-}
-
-Vec3 CompMultiply(const Vec3& v1, const Vec3& v2)
-{
-  Vec3 nrvo(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2]);
-  return nrvo;
-}
-
-Vec3 CompDivide(const Vec3& v1, const Vec3& v2)
-{
-  Vec3 nrvo(v1[0]/v2[0],v1[1]/v2[1],v1[2]/v2[2]);
-  return nrvo;
-}
-
-Vec3 operator*(const Vec3& v,const Mat3& m)
-{
-  Vec3 nrvo(v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0),
-      v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1),
-      v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2));
-  return nrvo;
-}
-
-Vec3 operator*(const Mat3& m, const Vec3& v)
-{
-  Vec3 nrvo(v[0]*m(0,0)+v[1]*m(0,1)+v[2]*m(0,2),
-      v[0]*m(1,0)+v[1]*m(1,1)+v[2]*m(1,2),
-      v[0]*m(2,0)+v[1]*m(2,1)+v[2]*m(2,2));
-  return nrvo;
-}
-
-Mat3 operator*(const Mat3& m1, const Mat3& m2)
-{
-  Mat3 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0)+m1(0,2)*m2(2,0),
-      m1(0,0)*m2(0,1)+m1(0,1)*m2(1,1)+m1(0,2)*m2(2,1),
-      m1(0,0)*m2(0,2)+m1(0,1)*m2(1,2)+m1(0,2)*m2(2,2),
-
-      m1(1,0)*m2(0,0)+m1(1,1)*m2(1,0)+m1(1,2)*m2(2,0),
-      m1(1,0)*m2(0,1)+m1(1,1)*m2(1,1)+m1(1,2)*m2(2,1),
-      m1(1,0)*m2(0,2)+m1(1,1)*m2(1,2)+m1(1,2)*m2(2,2),
-
-      m1(2,0)*m2(0,0)+m1(2,1)*m2(1,0)+m1(2,2)*m2(2,0),
-      m1(2,0)*m2(0,1)+m1(2,1)*m2(1,1)+m1(2,2)*m2(2,1),
-      m1(2,0)*m2(0,2)+m1(2,1)*m2(1,2)+m1(2,2)*m2(2,2));
-  return nrvo;
-}
-
 Real Comp(const Mat3& m,unsigned int i_in, unsigned int j_in)
 {
   return Minor(m,i_in,j_in)* static_cast<Real>( ((i_in+j_in)&0x1) ? -1 : 1);
@@ -264,25 +168,4 @@ Mat3 AxisRotation(const Vec3& axis, Real angle)
           xz-ca*xz-sa*y, yz-ca*yz+sa*x,zz+ca-ca*zz);
 }
 
-Vec3 Min(const Vec3& v1, const Vec3& v2)
-{
-  Vec3 nrvo(std::min(v1[0],v2[0]),
-            std::min(v1[1],v2[1]),
-            std::min(v1[2],v2[2]));
-  return nrvo;
-}
-
-Vec3 Max(const Vec3& v1, const Vec3& v2)
-{
-  Vec3 nrvo(std::max(v1[0],v2[0]),
-            std::max(v1[1],v2[1]),
-            std::max(v1[2],v2[2]));
-  return nrvo;
-}
-
-Real  Distance(const Vec3& p1, const Vec3& p2)
-{
-  return Length(p1-p2);
-}
-
 } // ns
diff --git a/modules/geom/src/vecmat3_op.hh b/modules/geom/src/vecmat3_op.hh
index 44c835768d3355e5c01e7785df971d9886162494..e4ca5a65b5eb0b2a9d11717fb68fa7a2c2c14a9d 100644
--- a/modules/geom/src/vecmat3_op.hh
+++ b/modules/geom/src/vecmat3_op.hh
@@ -23,48 +23,121 @@
 #include "constants.hh"
 
 #include <ost/geom/module_config.hh>
+#include <ost/geom/vec3.hh>
+#include <ost/geom/mat3.hh>
 
 namespace geom {
 
-class Vec3;
-class Mat3;
-
-//! returns length of vector
-Real DLLEXPORT_OST_GEOM Length(const Vec3& v);
 
 //! returns squared length of vector
-Real DLLEXPORT_OST_GEOM Length2(const Vec3& v);
+inline Real Length2(const Vec3& v)
+{
+  return v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
+}
 
-//! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Vec3& v1, const Vec3& v2, 
-                              Real ephilon=EPSILON);
+//! returns length of vector
+inline Real Length(const Vec3& v)
+{
+  return std::sqrt(Length2(v));
+}
+
+//! return true if components differ not more than epsilon
+inline bool Equal(const Vec3& v1, const Vec3& v2, 
+                  Real ephilon=EPSILON)
+{
+  return std::fabs(v1[0]-v2[0])<ephilon &&
+    std::fabs(v1[1]-v2[1])<ephilon &&
+    std::fabs(v1[2]-v2[2])<ephilon;
+}
 
 //! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Mat3& m1, const Mat3& m2, 
-                              Real ephilon=EPSILON);
+inline bool Equal(const Mat3& m1, const Mat3& m2, 
+                  Real ephilon=EPSILON)
+{
+  return std::fabs(m1(0,0)-m2(0,0))<ephilon &&
+    std::fabs(m1(0,1)-m2(0,1))<ephilon &&
+    std::fabs(m1(0,2)-m2(0,2))<ephilon &&
+    std::fabs(m1(1,0)-m2(1,0))<ephilon &&
+    std::fabs(m1(1,1)-m2(1,1))<ephilon &&
+    std::fabs(m1(1,2)-m2(1,2))<ephilon &&
+    std::fabs(m1(2,0)-m2(2,0))<ephilon &&
+    std::fabs(m1(2,1)-m2(2,1))<ephilon &&
+    std::fabs(m1(2,2)-m2(2,2))<ephilon;
+}
 
 //! vector dot product
-Real DLLEXPORT_OST_GEOM Dot(const Vec3& v1, const Vec3& v2);
+inline Real Dot(const Vec3& v1, const Vec3& v2)
+{
+  return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2];
+}
 
 //! Normalize
-Vec3 DLLEXPORT_OST_GEOM Normalize(const Vec3& v);
+inline Vec3 Normalize(const Vec3& v)
+{
+  Real l=Length(v);
+  if(l==0.0) {
+    return v;
+  }
+  return v/l;
+}
 
 //! vector cross product
-Vec3 DLLEXPORT_OST_GEOM Cross(const Vec3& v1, const Vec3& v2);
+inline Vec3 Cross(const Vec3& v1, const Vec3& v2)
+{
+  Vec3 nrvo(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 nrvo;
+}
 
 //! multiply each component of v1 with that of v2
-Vec3 DLLEXPORT_OST_GEOM CompMultiply(const Vec3& v1, const Vec3& v2);
+inline Vec3 CompMultiply(const Vec3& v1, const Vec3& v2)
+{
+  Vec3 nrvo(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2]);
+  return nrvo;
+}
 
 //! divide each component of v1 by that of v2
-Vec3 DLLEXPORT_OST_GEOM CompDivide(const Vec3& v1, const Vec3& v2);
+inline Vec3 CompDivide(const Vec3& v1, const Vec3& v2)
+{
+  Vec3 nrvo(v1[0]/v2[0],v1[1]/v2[1],v1[2]/v2[2]);
+  return nrvo;
+}
 
 //! vector matrix multiplication
-Vec3 DLLEXPORT_OST_GEOM operator*(const Vec3& v,const Mat3& m);
+inline Vec3 operator*(const Vec3& v,const Mat3& m)
+{
+  Vec3 nrvo(v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0),
+      v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1),
+      v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2));
+  return nrvo;
+}
+
 //! vector matrix multiplication
-Vec3 DLLEXPORT_OST_GEOM operator*(const Mat3& m, const Vec3& v);
+inline Vec3 operator*(const Mat3& m, const Vec3& v)
+{
+  Vec3 nrvo(v[0]*m(0,0)+v[1]*m(0,1)+v[2]*m(0,2),
+      v[0]*m(1,0)+v[1]*m(1,1)+v[2]*m(1,2),
+      v[0]*m(2,0)+v[1]*m(2,1)+v[2]*m(2,2));
+  return nrvo;
+}
 
 //! matrix matrix multiplication
-Mat3 DLLEXPORT_OST_GEOM operator*(const Mat3& m1, const Mat3& m2);
+inline Mat3 operator*(const Mat3& m1, const Mat3& m2)
+{
+  Mat3 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0)+m1(0,2)*m2(2,0),
+      m1(0,0)*m2(0,1)+m1(0,1)*m2(1,1)+m1(0,2)*m2(2,1),
+      m1(0,0)*m2(0,2)+m1(0,1)*m2(1,2)+m1(0,2)*m2(2,2),
+
+      m1(1,0)*m2(0,0)+m1(1,1)*m2(1,0)+m1(1,2)*m2(2,0),
+      m1(1,0)*m2(0,1)+m1(1,1)*m2(1,1)+m1(1,2)*m2(2,1),
+      m1(1,0)*m2(0,2)+m1(1,1)*m2(1,2)+m1(1,2)*m2(2,2),
+
+      m1(2,0)*m2(0,0)+m1(2,1)*m2(1,0)+m1(2,2)*m2(2,0),
+      m1(2,0)*m2(0,1)+m1(2,1)*m2(1,1)+m1(2,2)*m2(2,1),
+      m1(2,0)*m2(0,2)+m1(2,1)*m2(1,2)+m1(2,2)*m2(2,2));
+  return nrvo;
+}
 
 Mat3 DLLEXPORT_OST_GEOM Invert(const Mat3& m);
 Mat3 DLLEXPORT_OST_GEOM Transpose(const Mat3& m);
@@ -91,13 +164,28 @@ Mat3 DLLEXPORT_OST_GEOM AxisRotation(const Vec3& axis, Real angle);
 Vec3 DLLEXPORT_OST_GEOM OrthogonalVector(const Vec3& axis);
 
 //! returns std::min of each component
-Vec3 DLLEXPORT_OST_GEOM Min(const Vec3& v1, const Vec3& v2);
+inline Vec3 Min(const Vec3& v1, const Vec3& v2)
+{
+  Vec3 nrvo(std::min(v1[0],v2[0]),
+            std::min(v1[1],v2[1]),
+            std::min(v1[2],v2[2]));
+  return nrvo;
+}
 
 //! returns std::max of each component
-Vec3 DLLEXPORT_OST_GEOM Max(const Vec3& v1, const Vec3& v2);
+inline Vec3 Max(const Vec3& v1, const Vec3& v2)
+{
+  Vec3 nrvo(std::max(v1[0],v2[0]),
+            std::max(v1[1],v2[1]),
+            std::max(v1[2],v2[2]));
+  return nrvo;
+}
 
 //! returns the distance between two points
-Real DLLEXPORT_OST_GEOM Distance(const Vec3& p1, const Vec3& p2);
+inline Real Distance(const Vec3& p1, const Vec3& p2)
+{
+    return Length(p1-p2);
+}
 
 } // ns
 
diff --git a/modules/geom/src/vecmat4_op.cc b/modules/geom/src/vecmat4_op.cc
index a677c788229be3dd48834acd0ebd2e42dcc401bd..4567575413b9950f4e67d185018a6561b0fbaf56 100644
--- a/modules/geom/src/vecmat4_op.cc
+++ b/modules/geom/src/vecmat4_op.cc
@@ -23,93 +23,9 @@
 #include "vecmat3_op.hh"
 
 #include "exc.hh"
-#include "vec4.hh"
-#include "mat4.hh"
-#include "mat3.hh"
 
-namespace geom {
-
-Real Length(const Vec4& v)
-{
-  return std::sqrt(Length2(v));
-}
-
-Real Length2(const Vec4& v)
-{
-  return v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
-}
-
-bool Equal(const Vec4& v1, const Vec4& v2, Real epsilon)
-{
-  return std::fabs(v1[0]-v2[0])<epsilon &&
-    std::fabs(v1[1]-v2[1])<epsilon &&
-    std::fabs(v1[2]-v2[2])<epsilon &&
-    std::fabs(v1[3]-v2[3])<epsilon;
-}
-
-bool Equal(const Mat4& m1, const Mat4& m2, Real epsilon)
-{
-  return std::fabs(m1(0,0)-m2(0,0))<epsilon &&
-    std::fabs(m1(0,1)-m2(0,1))<epsilon &&
-    std::fabs(m1(0,2)-m2(0,2))<epsilon &&
-    std::fabs(m1(0,3)-m2(0,3))<epsilon &&
-    std::fabs(m1(1,0)-m2(1,0))<epsilon &&
-    std::fabs(m1(1,1)-m2(1,1))<epsilon &&
-    std::fabs(m1(1,2)-m2(1,2))<epsilon &&
-    std::fabs(m1(1,3)-m2(1,3))<epsilon &&
-    std::fabs(m1(2,0)-m2(2,0))<epsilon &&
-    std::fabs(m1(2,1)-m2(2,1))<epsilon &&
-    std::fabs(m1(2,2)-m2(2,2))<epsilon &&
-    std::fabs(m1(2,3)-m2(2,3))<epsilon &&
-    std::fabs(m1(3,0)-m2(3,0))<epsilon &&
-    std::fabs(m1(3,1)-m2(3,1))<epsilon &&
-    std::fabs(m1(3,2)-m2(3,2))<epsilon &&
-    std::fabs(m1(3,3)-m2(3,3))<epsilon;
-}
-
-Real Dot(const Vec4& v1, const Vec4& v2)
-{
-  return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3];
-}
-
-Vec4 Normalize(const Vec4& v)
-{
-  Real l=Length(v);
-  if(l==0.0) {
-    return v;
-  }
-  return v/l;
-}
-
-Vec4 CompMultiply(const Vec4& v1, const Vec4& v2)
-{
-  Vec4 nrvo(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2],v1[3]*v2[3]);
-  return nrvo;
-}
-
-Vec4 CompDivide(const Vec4& v1, const Vec4& v2)
-{
-  Vec4 nrvo(v1[0]/v2[0],v1[1]/v2[1],v1[2]/v2[2],v1[3]/v2[3]);
-  return nrvo;
-}
 
-Vec4 operator*(const Vec4& v,const Mat4& m)
-{
-  Vec4 nrvo(v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0)+v[3]*m(3,0),
-      v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1)+v[3]*m(3,1),
-      v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2)+v[3]*m(3,2),
-      v[0]*m(0,3)+v[1]*m(1,3)+v[2]*m(2,3)+v[3]*m(3,3));
-  return nrvo;
-}
-
-Vec4 operator*(const Mat4& m, const Vec4& v)
-{
-  Vec4 nrvo(v[0]*m(0,0)+v[1]*m(0,1)+v[2]*m(0,2)+v[3]*m(0,3),
-      v[0]*m(1,0)+v[1]*m(1,1)+v[2]*m(1,2)+v[3]*m(1,3),
-      v[0]*m(2,0)+v[1]*m(2,1)+v[2]*m(2,2)+v[3]*m(2,3),
-      v[0]*m(3,0)+v[1]*m(3,1)+v[2]*m(3,2)+v[3]*m(3,3));
-  return nrvo;
-}
+namespace geom {
 
 Real Comp(const Mat4& m, unsigned int i_in, unsigned int j_in)
 {
@@ -158,7 +74,6 @@ Mat4 Transpose(const Mat4& m)
   return r;
 }
 
-#if 1
 Mat4 Invert(const Mat4& in)
 {
   Mat4 m(in);
@@ -221,22 +136,6 @@ Mat4 Invert(const Mat4& in)
   return m;
 }
 
-#else
-
-Mat4 Invert(const Mat4& m)
-{
-  Mat4 r;
-  Real d=Det(m);
-  for (int i=0;i<4;i++){
-    for (int j=0;j<4;j++){
-      r(j,i)=Comp(m,i,j)/d;
-    }
-  }
-  return r;
-}
-
-#endif
-
 Mat4 operator*(const Mat4& m1, const Mat4& m2)
 {
   Mat4 nrvo(m1(0,0)*m2(0,0)+m1(0,1)*m2(1,0)+m1(0,2)*m2(2,0)+m1(0,3)*m2(3,0),
@@ -267,22 +166,4 @@ Real Angle(const Vec4& v1, const Vec4& v2)
   return std::acos(t);
 }
 
-Vec4 Min(const Vec4& v1, const Vec4& v2)
-{
-  Vec4 nrvo(std::min(v1[0],v2[0]),
-            std::min(v1[1],v2[1]),
-            std::min(v1[2],v2[2]),
-            std::min(v1[3],v2[3]));
-  return nrvo;
-}
-
-Vec4 Max(const Vec4& v1, const Vec4& v2)
-{
-  Vec4 nrvo(std::max(v1[0],v2[0]),
-            std::max(v1[1],v2[1]),
-            std::max(v1[2],v2[2]),
-            std::max(v1[3],v2[3]));
-  return nrvo;
-}
-
 } // ns
diff --git a/modules/geom/src/vecmat4_op.hh b/modules/geom/src/vecmat4_op.hh
index 501a09ebbe0f2684e383c2c5d8c675504b085629..7f7d6a2a95909da8dae244be343366106b7e5d8d 100644
--- a/modules/geom/src/vecmat4_op.hh
+++ b/modules/geom/src/vecmat4_op.hh
@@ -22,39 +22,104 @@
 #include "constants.hh"
 
 #include <ost/geom/module_config.hh>
+#include <ost/geom/vec4.hh>
+#include <ost/geom/vec3.hh>
+#include <ost/geom/mat4.hh>
+#include <ost/geom/mat3.hh>
 
 namespace geom {
 
-class Vec4;
-class Mat4;
-
-//! returns length of vector
-Real DLLEXPORT_OST_GEOM Length(const Vec4& v);
 
 //! returns squared length of vector
-Real DLLEXPORT_OST_GEOM Length2(const Vec4& v);
+inline Real Length2(const Vec4& v)
+{
+    return v[0]*v[0]+v[1]*v[1]+v[2]*v[2]+v[3]*v[3];
+}
+
+//! returns length of vector
+inline Real Length(const Vec4& v)
+{
+  return std::sqrt(Length2(v));
+}
 
 //! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Vec4& v1, const Vec4& v2, Real ephilon=EPSILON);
+inline bool Equal(const Vec4& v1, const Vec4& v2, Real epsilon=EPSILON)
+{
+  return std::fabs(v1[0]-v2[0])<epsilon &&
+    std::fabs(v1[1]-v2[1])<epsilon &&
+    std::fabs(v1[2]-v2[2])<epsilon &&
+    std::fabs(v1[3]-v2[3])<epsilon;
+}
 
 //! return true if components differ not more than ephilon
-bool DLLEXPORT_OST_GEOM Equal(const Mat4& m1, const Mat4& v2, Real ephilon=EPSILON);
+inline bool Equal(const Mat4& m1, const Mat4& m2, Real epsilon=EPSILON)
+{
+  return std::fabs(m1(0,0)-m2(0,0))<epsilon &&
+    std::fabs(m1(0,1)-m2(0,1))<epsilon &&
+    std::fabs(m1(0,2)-m2(0,2))<epsilon &&
+    std::fabs(m1(0,3)-m2(0,3))<epsilon &&
+    std::fabs(m1(1,0)-m2(1,0))<epsilon &&
+    std::fabs(m1(1,1)-m2(1,1))<epsilon &&
+    std::fabs(m1(1,2)-m2(1,2))<epsilon &&
+    std::fabs(m1(1,3)-m2(1,3))<epsilon &&
+    std::fabs(m1(2,0)-m2(2,0))<epsilon &&
+    std::fabs(m1(2,1)-m2(2,1))<epsilon &&
+    std::fabs(m1(2,2)-m2(2,2))<epsilon &&
+    std::fabs(m1(2,3)-m2(2,3))<epsilon &&
+    std::fabs(m1(3,0)-m2(3,0))<epsilon &&
+    std::fabs(m1(3,1)-m2(3,1))<epsilon &&
+    std::fabs(m1(3,2)-m2(3,2))<epsilon &&
+    std::fabs(m1(3,3)-m2(3,3))<epsilon;
+}
 
 //! vector dot product
-Real DLLEXPORT_OST_GEOM Dot(const Vec4& v1, const Vec4& v2);
-
-Vec4 DLLEXPORT_OST_GEOM Normalize(const Vec4& v);
+inline Real Dot(const Vec4& v1, const Vec4& v2)
+{
+    return v1[0]*v2[0]+v1[1]*v2[1]+v1[2]*v2[2]+v1[3]*v2[3];
+}
+
+inline Vec4 Normalize(const Vec4& v)
+{
+  Real l=Length(v);
+  if(l==0.0) {
+    return v;
+  }
+  return v/l;
+}
 
 //! multiply each component of v1 with that of v2
-Vec4 DLLEXPORT_OST_GEOM CompMultiply(const Vec4& v1, const Vec4& v2);
+inline Vec4 CompMultiply(const Vec4& v1, const Vec4& v2)
+{
+  Vec4 nrvo(v1[0]*v2[0],v1[1]*v2[1],v1[2]*v2[2],v1[3]*v2[3]);
+  return nrvo;
+}
 
 //! divide each component of v1 by that of v2
-Vec4 DLLEXPORT_OST_GEOM CompDivide(const Vec4& v1, const Vec4& v2);
+inline Vec4 CompDivide(const Vec4& v1, const Vec4& v2)
+{
+  Vec4 nrvo(v1[0]/v2[0],v1[1]/v2[1],v1[2]/v2[2],v1[3]/v2[3]);
+  return nrvo;
+}
 
 //! vector matrix multiplication
-Vec4 DLLEXPORT_OST_GEOM operator*(const Vec4& v,const Mat4& m);
+inline Vec4 operator*(const Vec4& v,const Mat4& m)
+{
+  Vec4 nrvo(v[0]*m(0,0)+v[1]*m(1,0)+v[2]*m(2,0)+v[3]*m(3,0),
+      v[0]*m(0,1)+v[1]*m(1,1)+v[2]*m(2,1)+v[3]*m(3,1),
+      v[0]*m(0,2)+v[1]*m(1,2)+v[2]*m(2,2)+v[3]*m(3,2),
+      v[0]*m(0,3)+v[1]*m(1,3)+v[2]*m(2,3)+v[3]*m(3,3));
+  return nrvo;
+}
+
 //! vector matrix multiplication
-Vec4 DLLEXPORT_OST_GEOM operator*(const Mat4& m, const Vec4& v);
+inline Vec4 operator*(const Mat4& m, const Vec4& v)
+{
+  Vec4 nrvo(v[0]*m(0,0)+v[1]*m(0,1)+v[2]*m(0,2)+v[3]*m(0,3),
+      v[0]*m(1,0)+v[1]*m(1,1)+v[2]*m(1,2)+v[3]*m(1,3),
+      v[0]*m(2,0)+v[1]*m(2,1)+v[2]*m(2,2)+v[3]*m(2,3),
+      v[0]*m(3,0)+v[1]*m(3,1)+v[2]*m(3,2)+v[3]*m(3,3));
+  return nrvo;
+}
 
 // algebraic complement
 Real DLLEXPORT_OST_GEOM Comp(const Mat4& m, unsigned int i, unsigned int j);
@@ -69,10 +134,24 @@ Mat4 DLLEXPORT_OST_GEOM operator*(const Mat4& m1, const Mat4& m2);
 Real DLLEXPORT_OST_GEOM Angle(const Vec4& v1, const Vec4& v2);
 
 //! returns std::min of each component
-Vec4 DLLEXPORT_OST_GEOM Min(const Vec4& v1, const Vec4& v2);
+inline Vec4 Min(const Vec4& v1, const Vec4& v2)
+{
+  Vec4 nrvo(std::min(v1[0],v2[0]),
+            std::min(v1[1],v2[1]),
+            std::min(v1[2],v2[2]),
+            std::min(v1[3],v2[3]));
+  return nrvo;
+}
 
 //! returns std::max of each component
-Vec4 DLLEXPORT_OST_GEOM Max(const Vec4& v1, const Vec4& v2);
+inline Vec4  Max(const Vec4& v1, const Vec4& v2)
+{
+  Vec4 nrvo(std::max(v1[0],v2[0]),
+            std::max(v1[1],v2[1]),
+            std::max(v1[2],v2[2]),
+            std::max(v1[3],v2[3]));
+  return nrvo;
+}
 
 } // ns
 
diff --git a/modules/gfx/src/impl/simple_renderer.cc b/modules/gfx/src/impl/simple_renderer.cc
index 4e228f0dc06f9e7fee003411bbf02f241de77e26..bfb30d3133648b5216a4916c006a17daea51ee56 100644
--- a/modules/gfx/src/impl/simple_renderer.cc
+++ b/modules/gfx/src/impl/simple_renderer.cc
@@ -99,7 +99,7 @@ geom::Vec3 SimpleRenderer::GetDoubleBondPlane(mol::BondHandle b)
   GetBondPartnerNormal(vec, n, bond_vec, partner2, a1, a2);
   if(n==0) {
     //set vec into x direction if no partners are present
-    vec.SetX(1);
+    vec.x=1;
   }
   vec=geom::Normalize(geom::Cross(vec,bond_vec));
   return vec;
@@ -159,7 +159,7 @@ void SimpleRenderer::PrepareRendering(GfxView& view, IndexedVertexArray& va)
         if(it->bond.GetBondOrder()==2) {
           vec=GetDoubleBondPlane(it->bond);
         } else {
-          vec.SetX(1);
+          vec.x=1;
         }
         Real dbd=options_->GetBondOrderDistance();
         VertexID id00 = va.Add(p0+dbd*vec,geom::Vec3(),it->atom1->color);
diff --git a/modules/gui/src/scene_selection.cc b/modules/gui/src/scene_selection.cc
index 72e8f4f55b0bf66f36c7c9caaf9c31bd7f7eee31..b408907f4282f62bf1600abb06b6d1a8305d93c3 100644
--- a/modules/gui/src/scene_selection.cc
+++ b/modules/gui/src/scene_selection.cc
@@ -110,9 +110,9 @@ void SceneSelection::CenterOnObjects() {
     }
   }
   if(changed){
-    gfx::Scene::Instance().SetCenter(geom::Vec3(((max.GetX()+min.GetX())/2),
-                                                ((max.GetY()+min.GetY())/2),
-                                                ((max.GetZ()+min.GetZ())/2)));
+    gfx::Scene::Instance().SetCenter(geom::Vec3(((max.x+min.x)/2),
+                                                ((max.y+min.y)/2),
+                                                ((max.z+min.z)/2)));
   }
 }
 
diff --git a/modules/io/src/mol/dcd_io.cc b/modules/io/src/mol/dcd_io.cc
index 0ac5a745632802c16be0b7a49a0e27c3aba9bb1c..29cbb7bd280b93437f00c072d5efba011846a19e 100644
--- a/modules/io/src/mol/dcd_io.cc
+++ b/modules/io/src/mol/dcd_io.cc
@@ -163,7 +163,7 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
     if(swap_flag) swap_float(&xlist[0],xlist.size());
     if(gap_flag) ff.read(dummy,sizeof(dummy));
     for(uint j=0;j<clist.size();++j) {
-      clist[j].SetX(xlist[j]);
+      clist[j].x=xlist[j];
     }
 
     // y coord
@@ -172,7 +172,7 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
     if(swap_flag) swap_float(&xlist[0],xlist.size());
     if(gap_flag) ff.read(dummy,sizeof(dummy));
     for(uint j=0;j<clist.size();++j) {
-      clist[j].SetY(xlist[j]);
+      clist[j].y=xlist[j];
     }
 
     // z coord
@@ -181,7 +181,7 @@ mol::CoordGroupHandle load_dcd(const mol::AtomHandleList& alist2,
     if(swap_flag) swap_float(&xlist[0],xlist.size());
     if(gap_flag) ff.read(dummy,sizeof(dummy));
     for(uint j=0;j<clist.size();++j) {
-      clist[j].SetZ(xlist[j]);
+      clist[j].z=xlist[j];
     }
 
     cg.AddFrame(clist);
diff --git a/modules/io/src/mol/entity_io_crd_handler.cc b/modules/io/src/mol/entity_io_crd_handler.cc
index 2fd39293ff0fe11c929e97915a59735edb034742..1b093be94b9478279763131256b9eb29f2b28a6b 100644
--- a/modules/io/src/mol/entity_io_crd_handler.cc
+++ b/modules/io/src/mol/entity_io_crd_handler.cc
@@ -186,9 +186,9 @@ bool CRDWriter::VisitAtom(const mol::AtomHandle& atom)
               << format("%5i") % res.GetNumber() << " "
               << format("%4s") % res.GetKey() << " "
               << format("%-4s") % atom.GetName()
-              << format("%10.5f") % atom.GetPos().GetX()
-              << format("%10.5f") % atom.GetPos().GetY()
-              << format("%10.5f") % atom.GetPos().GetZ() << " "
+              << format("%10.5f") % atom.GetPos().x
+              << format("%10.5f") % atom.GetPos().y
+              << format("%10.5f") % atom.GetPos().z << " "
               << format("%-4s") % e_name << " "
               << format("%-5i") % res.GetNumber() << " "
               << format("%8.5f") % atom.GetBFactor()