diff --git a/t.cc b/t.cc
index 0ae16fb..789de85 100644
--- a/t.cc
+++ b/t.cc
@@ -3715,6 +3715,16 @@ cout <<t;
 fourindex_dense<nosymmetry,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+const fourindex_dense<nosymmetry,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i) 
+for(int j=1; j<=n; ++j) 
+for(int k=1; k<=n; ++k) 
+for(int l=1; l<=n; ++l) 
+	{
+	cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+	if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+	}
 }
 
 if(0)
@@ -3728,6 +3738,18 @@ cout <<t;
 fourindex_dense<twoelectronrealmullikanAB,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<twoelectronrealmullikanAB,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n; ++k)
+for(int l=1; l<=n; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+        if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+        if(fc(i,j,k,l) != tc(l,k,j,i) ) laerror("operator() mismatch");
+        }
 }
 
 if(0)
@@ -3741,6 +3763,18 @@ cout <<t;
 fourindex_dense<twoelectronrealmullikan,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<twoelectronrealmullikan,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n; ++k)
+for(int l=1; l<=n; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<l<<" "<<k<<" "<<j<<" "<<i<<" : "<<tc(l,k,j,i)<<endl;
+        if(fc(i,j,k,l) != tc(l,k,j,i) ) laerror("operator() mismatch");
+        if(fc(i,j,k,l) != tc(i,j,k,l) ) laerror("operator() mismatch");
+        }
 }
 
 if(0)
@@ -3754,6 +3788,18 @@ cout <<t;
 fourindex_dense<T2IjAb_aces,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<T2IjAb_aces,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n+1; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n+3; ++k)
+for(int l=1; l<=n+2; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+        if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+        }
+
 }
 
 if(0)
@@ -3767,6 +3813,18 @@ cout <<t;
 fourindex_dense<T2ijab_aces,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<T2ijab_aces,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n+2; ++k)
+for(int l=1; l<=n+2; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+        if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+        }
+
 }
 
 if(0)
@@ -3780,6 +3838,18 @@ cout <<t;
 fourindex_dense<antisymtwoelectronrealdiracAB,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<antisymtwoelectronrealdiracAB,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n; ++k)
+for(int l=1; l<=n; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+       if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+       if(fc(i,j,k,l) != tc(i,j,k,l) ) laerror("operator() mismatch");
+        }
 }
 
 
@@ -3794,11 +3864,23 @@ cout <<t;
 fourindex_dense<antisymtwoelectronrealdirac,double,int> ff;
 tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
+
+const fourindex_dense<antisymtwoelectronrealdirac,double,int> fc(f);
+const Tensor<double> tc(t);
+for(int i=1; i<=n; ++i)
+for(int j=1; j<=n; ++j)
+for(int k=1; k<=n; ++k)
+for(int l=1; l<=n; ++l)
+        {
+        cout<< "4ind "<<i<<" "<<j<<" "<<k<<" "<<l<<" : " <<fc(i,j,k,l) <<" = tensor " <<k<<" "<<l<<" "<<i<<" "<<j<<" : "<<tc(k,l,i,j)<<endl;
+       if(fc(i,j,k,l) != tc(k,l,i,j) ) laerror("operator() mismatch");
+       if(fc(i,j,k,l) != tc(i,j,k,l) ) laerror("operator() mismatch");
+        }
+
 }
 
 
 
-//@@@check operator() for all fourindex-tensor conversions
 
 
 }//main
diff --git a/tensor.h b/tensor.h
index 6683f3c..ce1590a 100644
--- a/tensor.h
+++ b/tensor.h
@@ -275,6 +275,8 @@ void tensor2fourindex(const Tensor<T> &t, fourindex_dense<S,T,I> &f);
 
 //conversions from/to fourindex specialized by symmetry type
 //NOTE also different index order due to tensor's leftmost index being the least significant
+//
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<nosymmetry,T,I> &f)
 {
@@ -307,6 +309,7 @@ f=fourindex_dense<nosymmetry,T,I>(range,NRMat<T>(t.data,range*range,range*range)
 }
 
 
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j)=(l,k,j,i) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<twoelectronrealmullikanAB,T,I> &f)
 {
@@ -340,6 +343,7 @@ f=fourindex_dense<twoelectronrealmullikanAB,T,I>(NRMat<T>(t.data,range*(range+1)
 }
 
 
+//operator() : (i,j,k,l) in fourindex -> (i,j,k,l)=(l,k,j,i) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<twoelectronrealmullikan,T,I> &f)
 {
@@ -375,6 +379,7 @@ f=fourindex_dense<twoelectronrealmullikan,T,I>(NRSMat<T>(mat)); //symmetrize mat
 }
 
 
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<T2IjAb_aces,T,I> &f)
 {
@@ -385,10 +390,10 @@ for(int i=0; i<4; ++i)
 	shape[i].symmetry=0;
 	shape[i].offset=1;
 	}
-shape[3].range=f.noca;
-shape[2].range=f.nocb;
-shape[1].range=f.nvra;
-shape[0].range=f.nvrb;
+shape[2].range=f.noca;
+shape[3].range=f.nocb;
+shape[0].range=f.nvra;
+shape[1].range=f.nvrb;
 NRVec<T> data(f);
 return Tensor<T>(shape,data);
 }
@@ -397,10 +402,10 @@ template<typename T, typename I>
 void tensor2fourindex(const Tensor<T> &t, fourindex_dense<T2IjAb_aces,T,I> &f)
 {
 if(t.rank()!=4 ||t.shape.size()!=4) laerror("wrong rank/shape in tensor2fourindex");
-int noca = t.shape[3].range;
-int nocb = t.shape[2].range;
-int nvra = t.shape[1].range;
-int nvrb = t.shape[0].range;
+int noca = t.shape[2].range;
+int nocb = t.shape[3].range;
+int nvra = t.shape[0].range;
+int nvrb = t.shape[1].range;
 
 int offset=t.shape[0].offset;
 for(int i=0; i<t.shape.size(); ++i)
@@ -413,7 +418,7 @@ f=fourindex_dense<T2IjAb_aces,T,I>(noca,nocb,nvra,nvrb,mat);
 }
 
 
-
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<T2ijab_aces,T,I> &f)
 {
@@ -424,8 +429,8 @@ for(int i=0; i<2; ++i)
 	shape[i].symmetry= -1;
 	shape[i].offset=1;
 	}
-shape[1].range=f.nocc;
 shape[0].range=f.nvrt;
+shape[1].range=f.nocc;
 NRVec<T> data(f);
 return Tensor<T>(shape,data);
 }
@@ -448,7 +453,7 @@ f=fourindex_dense<T2ijab_aces,T,I>(nocc,nvrt,mat);
 }
 
 
-
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j)=(i,j,k,l) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<antisymtwoelectronrealdiracAB,T,I> &f)
 {
@@ -483,7 +488,7 @@ f=fourindex_dense<antisymtwoelectronrealdiracAB,T,I>(range,NRSMat<T>(mat));
 }
 
 
-
+//operator() : (i,j,k,l) in fourindex -> (k,l,i,j)=(i,j,k,l) in tensor
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<antisymtwoelectronrealdirac,T,I> &f)
 {