diff --git a/examples/qm/2d-coulomb-small-TEN-2.py b/examples/qm/2d-coulomb-small-TEN-2.py
index c94052bdb..a360fc7e5 100644
--- a/examples/qm/2d-coulomb-small-TEN-2.py
+++ b/examples/qm/2d-coulomb-small-TEN-2.py
@@ -35,7 +35,7 @@
 		[Ip2_QMParameters_QMParametersCoulomb_QMIPhysCoulomb()],
 		[Law2_QMIGeom_QMIPhysCoulomb()]
 	),
-	SchrodingerKosloffPropagator(steps=-1,threadNum=4),
+	SchrodingerKosloffPropagator(steps=-1,threadNum=8),
 ]
 
 partsScale = 20000
@@ -90,24 +90,24 @@
 #O.save('/tmp/a.xml.bz2');
 #o.run(100000); o.wait(); print o.iter/o.realtime,'iterations/sec'
 
-try:
-	from yade import qt
-	qt.Controller()
-	qt.controller.setViewAxes(dir=(0,1,0),up=(0,0,1))
-	qt.controller.setWindowTitle(sys.argv[0])
-	qt.Renderer().blinkHighlight=False
-	qt.View()
-        qt.Renderer().light2Pos=[Pot_x,Pot_y,30]
-	qt.views()[0].center(False,250) # median=False, suggestedRadius = 5
+## try:
+## 	from yade import qt
+## 	qt.Controller()
+## 	qt.controller.setViewAxes(dir=(0,1,0),up=(0,0,1))
+## 	qt.controller.setWindowTitle(sys.argv[0])
+## 	qt.Renderer().blinkHighlight=False
+## 	qt.View()
+##         qt.Renderer().light2Pos=[Pot_x,Pot_y,30]
+## 	qt.views()[0].center(False,250) # median=False, suggestedRadius = 5
+## 
+## except ImportError:
+## 	pass
 
-except ImportError:
-	pass
 
+O.dt=0.0000001
+for i in range(51):
+	O.step()
+	O.dt=200
+	if(i%5==0):
+		O.save(str(sys.argv[0])+"_"+str(O.iter)+".yade.gz")
 
-#O.dt=0.0000001
-#for i in range(51):
-#	O.step()
-#	O.dt=200
-#	if(i%5==0):
-#		O.save(str(sys.argv[0])+"_"+str(O.iter)+".yade.bz2")
-#
diff --git a/examples/qm/2d-hydrogeneigen+positron-small-TEN-2small.py b/examples/qm/2d-hydrogeneigen+positron-small-TEN-2small.py
index 6e2399a03..157646395 100644
--- a/examples/qm/2d-hydrogeneigen+positron-small-TEN-2small.py
+++ b/examples/qm/2d-hydrogeneigen+positron-small-TEN-2small.py
@@ -46,7 +46,7 @@
                  ,Ip2_QMParticleCoulomb_QMParametersCoulomb_QMIPhysCoulombParticleInPotential()],
 		[Law2_QMIGeom_QMIPhysCoulombParticles(),Law2_QMIGeom_QMIPhysCoulombParticleInPotential()]
 	),
-	SchrodingerKosloffPropagator(FIXMEatomowe_MASS=1.0,steps=-1,virialCheck=False,printIter=20,doCopyTable=False,threadNum=8),
+	SchrodingerKosloffPropagator(FIXMEatomowe_MASS=1.0,steps=-1,virialCheck=False,printIter=1,doCopyTable=False,threadNum=16),
 	SchrodingerAnalyticPropagator(),
 ]
 
@@ -159,9 +159,9 @@ def placeDraw(di,i,j):
 	qt.Renderer().blinkHighlight=False
 	qt.Renderer().light1Pos=Vector3( 1175,1130,500)
 	qt.Renderer().light2Pos=Vector3(-1130, 575,230)
-	qt.View()
+	#qt.View()
 	#qt.Renderer().light2Pos=Vector3(Pot_x,Pot_y,30)
-	qt.views()[0].center(False,size1d*1.5) # median=False, suggestedRadius = 5
+	#qt.views()[0].center(False,size1d*1.5) # median=False, suggestedRadius = 5
 
 except ImportError:
 	pass
diff --git a/examples/qm/3d-coulomb-eigen-wavefunction.py b/examples/qm/3d-coulomb-eigen-wavefunction.py
index ff958ecbe..c97dda537 100644
--- a/examples/qm/3d-coulomb-eigen-wavefunction.py
+++ b/examples/qm/3d-coulomb-eigen-wavefunction.py
@@ -34,7 +34,7 @@
 		[Ip2_QMParameters_QMParametersCoulomb_QMIPhysCoulomb()],
 		[Law2_QMIGeom_QMIPhysCoulomb()]
 	),
-	SchrodingerKosloffPropagator(),
+	SchrodingerKosloffPropagator(threadNum=8),
 	SchrodingerAnalyticPropagator(),
 ]
 
diff --git a/lib/base/NDimTable.hpp b/lib/base/NDimTable.hpp
index fecc14153..970bee8b7 100644
--- a/lib/base/NDimTable.hpp
+++ b/lib/base/NDimTable.hpp
@@ -25,6 +25,7 @@
 #include <iomanip>
 #include <boost/lexical_cast.hpp>
 #include <Eigen/Core>
+#include <parallel/algorithm>
 
 #include <boost/thread/mutex.hpp>
 #include <boost/serialization/nvp.hpp>
@@ -343,18 +344,18 @@ std::cerr << "destructor                          : " << --ZZ::NDimTable_Instanc
 		// elementwise operations
 		NDimTable& operator  = (const NDimTable& )=default;
 		NDimTable& operator  = (      NDimTable&&)=default;
-		NDimTable& operator -  (          ) {std::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return  -v;});  return *this;}; 
-		NDimTable& operator += (const K& k) {std::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v+k;});  return *this;};
-		NDimTable& operator -= (const K& k) {std::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v-k;});  return *this;}; 
-		NDimTable& operator *= (const K& k) {std::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v*k;});  return *this;}; 
-		NDimTable& operator /= (const K& k) {std::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v/k;});  return *this;}; 
-
-		template<typename L> NDimTable& operator += (const NDimTable<L>& T) {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v+l;});return *this;};
-		template<typename L> NDimTable& operator -= (const NDimTable<L>& T) {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v-l;});return *this;};
+		NDimTable& operator -  (          ) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return  -v;});  return *this;};
+		NDimTable& operator += (const K& k) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v+k;});  return *this;};
+		NDimTable& operator -= (const K& k) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v-k;});  return *this;};
+		NDimTable& operator *= (const K& k) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v*k;});  return *this;};
+		NDimTable& operator /= (const K& k) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[k](K& v){return v/k;});  return *this;};
+
+		template<typename L> NDimTable& operator += (const NDimTable<L>& T) {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v+l;});return *this;};
+		template<typename L> NDimTable& operator -= (const NDimTable<L>& T) {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v-l;});return *this;};
 		// !!!!!!!!!!!
 		// !IMPORTANT! operators *= and /= implement  http://en.wikipedia.org/wiki/Hadamard_product_%28matrices%29
-		template<typename L> NDimTable& operator *= (const NDimTable<L>& T) {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v*l;});return *this;};
-		template<typename L> NDimTable& operator /= (const NDimTable<L>& T) {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v/l;});return *this;}; 
+		template<typename L> NDimTable& operator *= (const NDimTable<L>& T) {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v*l;});return *this;};
+		template<typename L> NDimTable& operator /= (const NDimTable<L>& T) {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[](K& v,const L& l){return v/l;});return *this;};
 
 		// FIXME: should be 'K'-type not 'double'-type. But min(), max() works only with real numbers. 
 		// FIXME: this is because potential should be real valued (but isn't)
@@ -363,21 +364,21 @@ std::cerr << "destructor                          : " << --ZZ::NDimTable_Instanc
 		not_complex maxReal() const {not_complex ret(std::real(this->front())); for(K v : (*this)){ret  = std::max(std::real(v),ret);}; return ret;};
 		// !!!!!!!!!!!
 		// !IMPORTANT! for effciency, these do not copy construct new data, they modify in-place!
-		NDimTable& abs()           {std::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::abs(v    );}); return *this;};
-		NDimTable& pow(const K& k) {std::transform(this->begin(),this->end(),this->begin(),[k](K& v){return std::pow(v,k  );}); return *this;};
-		NDimTable& sqrt()          {std::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::sqrt(v   );}); return *this;};
-		NDimTable& conj()          {std::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::conj(v   );}); return *this;};
+		NDimTable& abs()           {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::abs(v    );}); return *this;};
+		NDimTable& pow(const K& k) {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[k](K& v){return std::pow(v,k  );}); return *this;};
+		NDimTable& sqrt()          {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::sqrt(v   );}); return *this;};
+		NDimTable& conj()          {__gnu_parallel::transform(this->begin(),this->end(),this->begin(),[ ](K& v){return std::conj(v   );}); return *this;};
 
 		template<typename L> NDimTable& mult2Add(const NDimTable<L>& T,const K& k)
-		                           {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v+l*k;});return *this;};
+		                           {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v+l*k;});return *this;};
 		template<typename L> NDimTable& mult2Sub(const NDimTable<L>& T,const K& k)
-		                           {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v-l*k;});return *this;};
+		                           {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v-l*k;});return *this;};
 		template<typename L> NDimTable& mult1Sub(const K& k,const NDimTable<L>& T)
-		                           {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v*k-l;});return *this;};
+		                           {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v*k-l;});return *this;};
 		template<typename L> NDimTable& multMult(const NDimTable<L>& T,const K& k)
-		                           {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v*l*k;});return *this;};
+		                           {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[k](K& v,const L& l){return v*l*k;});return *this;};
 		template<typename L> NDimTable& mult1Mult2Add(const K& k1,const NDimTable<L>& T,const K& k2)
-		                           {std::transform(this->begin(),this->end(),T.begin(),this->begin(),[k1,k2](K& v,const L& l){return v*k1+l*k2;});return *this;};
+		                           {__gnu_parallel::transform(this->begin(),this->end(),T.begin(),this->begin(),[k1,k2](K& v,const L& l){return v*k1+l*k2;});return *this;};
 
 		// // contractions (returns new container of different dimension, or works on a provided container of expected dimension)
 		// //