diff --git a/t.cc b/t.cc
index b27c65a..dd1612c 100644
--- a/t.cc
+++ b/t.cc
@@ -4597,6 +4597,8 @@ Tensor<double> x(shape); x.randomize(1.);
 x.defaultnames();
 cout <<"x= "<<x.shape << " "<<x.names<<endl;
 
+cout <<"indexmatrix of x = "<<x.indexmatrix();
+
 NRVec<INDEXGROUP> yshape(2);
         yshape[0].number=1; 
         yshape[0].symmetry=0;
diff --git a/tensor.cc b/tensor.cc
index c2704df..d2115dc 100644
--- a/tensor.cc
+++ b/tensor.cc
@@ -417,7 +417,7 @@ calcsize();
 
 
 template<typename T>
-void loopingroups(Tensor<T> &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *))
+void loopingroups(Tensor<T> &t, int ngroup, int igroup, T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, T *), bool skipzeros)
 {
 LA_index istart,iend;
 const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[ngroup];
@@ -443,6 +443,8 @@ switch(sh->symmetry)
 
 for(LA_index i = istart; i<=iend; ++i)
 	{
+	if(skipzeros && i==0) continue;
+
 	I[ngroup][igroup]=i;
 	if(ngroup==0 && igroup==0) 
 		{
@@ -460,14 +462,14 @@ for(LA_index i = istart; i<=iend; ++i)
 			const INDEXGROUP *sh2 = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[newngroup];
 			newigroup=sh2->number-1;
 			}
-		loopingroups(t,newngroup,newigroup,p,I,callback);
+		loopingroups(t,newngroup,newigroup,p,I,callback,skipzeros);
 		}
 	}
 }
 
 
 template<typename T>
-void Tensor<T>::loopover(void (*callback)(const SUPERINDEX &, T *))
+void Tensor<T>::loopover(void (*callback)(const SUPERINDEX &, T *), bool skipzeros)
 {
 SUPERINDEX I(shape.size());
 for(int i=0; i<I.size(); ++i)
@@ -479,12 +481,12 @@ for(int i=0; i<I.size(); ++i)
 T *pp=&data[0];
 int ss=shape.size()-1;
 const INDEXGROUP *sh = &(* const_cast<const NRVec<INDEXGROUP> *>(&shape))[ss];
-loopingroups(*this,ss,sh->number-1,&pp,I,callback);
+loopingroups(*this,ss,sh->number-1,&pp,I,callback,skipzeros);
 }
 
 
 template<typename T>
-void constloopingroups(const Tensor<T> &t, int ngroup, int igroup, const T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, const T *))
+void constloopingroups(const Tensor<T> &t, int ngroup, int igroup, const T **p, SUPERINDEX &I, void (*callback)(const SUPERINDEX &, const T *), bool skipzeros)
 {
 LA_index istart,iend;
 const INDEXGROUP *sh = &t.shape[ngroup];
@@ -510,6 +512,8 @@ switch(sh->symmetry)
 
 for(LA_index i = istart; i<=iend; ++i)
 	{
+        if(skipzeros && i==0) continue;
+
 	I[ngroup][igroup]=i;
 	if(ngroup==0 && igroup==0) 
 		{
@@ -527,14 +531,14 @@ for(LA_index i = istart; i<=iend; ++i)
 			const INDEXGROUP *sh2 = &(* const_cast<const NRVec<INDEXGROUP> *>(&t.shape))[newngroup];
 			newigroup=sh2->number-1;
 			}
-		constloopingroups(t,newngroup,newigroup,p,I,callback);
+		constloopingroups(t,newngroup,newigroup,p,I,callback,skipzeros);
 		}
 	}
 }
 
 
 template<typename T>
-void Tensor<T>::constloopover(void (*callback)(const SUPERINDEX &, const T *)) const
+void Tensor<T>::constloopover(void (*callback)(const SUPERINDEX &, const T *), bool skipzeros) const
 {
 SUPERINDEX I(shape.size());
 for(int i=0; i<I.size(); ++i)
@@ -546,7 +550,28 @@ for(int i=0; i<I.size(); ++i)
 const T *pp=&data[0];
 int ss=shape.size()-1;
 const INDEXGROUP *sh = &shape[ss];
-constloopingroups(*this,ss,sh->number-1,&pp,I,callback);
+constloopingroups(*this,ss,sh->number-1,&pp,I,callback,skipzeros);
+}
+
+static INDEXMATRIX *indexmat_p;
+static LA_largeindex indexmat_row;
+
+template<typename T>
+static void indexmatrix_callback(const SUPERINDEX &I, const T *p)
+{
+FLATINDEX f=superindex2flat(I);
+indexmat_p->rowset(f,indexmat_row++);
+}
+
+template<typename T>
+INDEXMATRIX  Tensor<T>::indexmatrix() const
+{
+INDEXMATRIX r;
+r.resize(size(),rank());
+indexmat_p = &r;
+indexmat_row = 0;
+constloopover(indexmatrix_callback<T>,false);
+return r;
 }
 
 
@@ -2439,6 +2464,27 @@ return true;
 }
 
 
+bool zero_in_index(const FLATINDEX &I)
+{
+for(int i=0; i<I.size(); ++i) if(I[i]==0) return true;
+return false;
+}
+
+bool zero_in_index(const INDEXMATRIX &m, const LA_largeindex row) 
+{
+const LA_index *p = &m(row,0);
+for(int i=0; i<m.ncols(); ++i) if(p[i]==0) return true;
+return false;
+}
+
+
+
+bool zero_in_index(const SUPERINDEX &I)
+{
+for(int i=0; i<I.size(); ++i) if(zero_in_index(I[i])) return true;
+return false;
+}
+
 
 template class Tensor<double>;
 template class Tensor<std::complex<double> >;
diff --git a/tensor.h b/tensor.h
index 3df491f..0a1f13e 100644
--- a/tensor.h
+++ b/tensor.h
@@ -190,8 +190,9 @@ class LA_traits<INDEXGROUP> {
 
 
 typedef NRVec<LA_index> FLATINDEX; //all indices but in a single vector
-typedef NRVec<NRVec<LA_index> > SUPERINDEX; //all indices in the INDEXGROUP structure
+typedef NRVec<FLATINDEX> SUPERINDEX; //all indices in the INDEXGROUP structure
 typedef NRVec<LA_largeindex> GROUPINDEX; //set of indices in the symmetry groups
+typedef NRMat<LA_index> INDEXMATRIX; //list of FLATINDEXes (rows of the matrix) of all tensor elements - convenient to be able to run over the whole tensor in a for loop rather than via recursive loopovers with a callback
 struct  INDEX
 {
 int group;
@@ -230,6 +231,11 @@ int flatposition(int group, int index, const NRVec<INDEXGROUP> &shape);
 int flatposition(const INDEX &i, const NRVec<INDEXGROUP> &shape); //position of that index in FLATINDEX
 INDEX indexposition(int flatindex, const NRVec<INDEXGROUP> &shape); //inverse to flatposition
 
+//useful for negative offsets and 0 index excluded
+bool zero_in_index(const FLATINDEX &);
+bool zero_in_index(const SUPERINDEX &);
+bool zero_in_index(const INDEXMATRIX &, const LA_largeindex row);
+
 FLATINDEX superindex2flat(const SUPERINDEX &I);
 
 template<typename T>
@@ -361,11 +367,14 @@ public:
 	bool fulfills_hermiticity() const; //check it is so
 	inline void randomize(const typename LA_traits<T>::normtype &x) {data.randomize(x); enforce_hermiticity();};
 
-	void loopover(void (*callback)(const SUPERINDEX &, T *)); //loop over all elements
-	void constloopover(void (*callback)(const SUPERINDEX &, const T *)) const; //loop over all elements
+
+	void loopover(void (*callback)(const SUPERINDEX &, T *), bool skipzeros=false); //loop over all elements, optionally skip zero indices (i.e. run over ...-2,-1,1,2...) which is useful for special applications
+	void constloopover(void (*callback)(const SUPERINDEX &, const T *), bool skipzeros=false) const; //loop over all elements
 	void grouploopover(void (*callback)(const GROUPINDEX &, T *)); //loop over all elements disregarding the internal structure of index groups
 	void constgrouploopover(void (*callback)(const GROUPINDEX &, const T *)) const; //loop over all elements disregarding the internal structure of index groups
 
+	INDEXMATRIX indexmatrix() const; //get indexmatrix - rows store FLATINDEXes matching data[]
+
 	Tensor permute_index_groups(const NRPerm<int> &p) const; //rearrange the tensor storage permuting index groups as a whole: result_i = source_p_i
 	Tensor permute_index_groups(const NRVec<INDEXNAME> &names) const; //permute to requested order of group's first indices (or permute individual indices of a flat tensor)