tensor-fourindex conversions for some more symemtries

fourindex.h

@@ -1182,12 +1182,13 @@ void print(std::ostream &out) const
 template<class T, class I>
 class fourindex_dense<T2IjAb_aces,T,I> : public NRMat<T> {
 protected:
-	unsigned int noca,nocb,nvra,nvrb;
 friend class explicit_t2;
 public:
+	unsigned int noca,nocb,nvra,nvrb;
         fourindex_dense(): NRMat<T>() {noca=nocb=nvra=nvrb=0;};
 	void resize(const int nocca, const int noccb, const int nvrta, const int nvrtb) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb; (*this).NRMat<T>::resize(nocca*noccb,nvrta*nvrtb);};
         explicit fourindex_dense(const int nocca, const int noccb, const int nvrta, const int nvrtb): NRMat<T>(nocca*noccb,nvrta*nvrtb) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb;};
+        explicit fourindex_dense(const int nocca, const int noccb, const int nvrta, const int nvrtb, const NRMat<T> &data): NRMat<T>(data) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb; if(data.nrows()!=nocca*noccb||data.ncols()!=nvrta*nvrtb) laerror("data size mismatch in fourindex_dense T2IjAb_aces constructor");};
 
 //here i,a are alpha j,b beta
         inline T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b)

t.cc

@@ -3704,7 +3704,7 @@ cout <<"Error = "<<(x1-y).norm()<<endl;
 }
 
 
-if(1)
+if(0)
 {
 int n;
 cin>>n;
@@ -3717,5 +3717,47 @@ tensor2fourindex(t,ff);
 cout <<"Error = "<<(f-ff).norm()<<endl;
 }
 
+if(0)
+{
+int n;
+cin>>n;
+fourindex_dense<twoelectronrealmullikanAB,double,int> f(n);
+f.randomize(1.);
+Tensor<double> t=fourindex2tensor(f);
+cout <<t;
+fourindex_dense<twoelectronrealmullikanAB,double,int> ff;
+tensor2fourindex(t,ff);
+cout <<"Error = "<<(f-ff).norm()<<endl;
+}
+
+if(0)
+{
+int n;
+cin>>n;
+fourindex_dense<twoelectronrealmullikan,double,int> f(n);
+f.randomize(1.);
+Tensor<double> t=fourindex2tensor(f);
+cout <<t;
+fourindex_dense<twoelectronrealmullikan,double,int> ff;
+tensor2fourindex(t,ff);
+cout <<"Error = "<<(f-ff).norm()<<endl;
+}
+
+if(1)
+{
+//check also index order
+int n;
+cin>>n;
+fourindex_dense<T2IjAb_aces,double,int> f(n+1,n,n+3,n+2);
+f.randomize(1.);
+Tensor<double> t=fourindex2tensor(f);
+cout <<t;
+fourindex_dense<T2IjAb_aces,double,int> ff;
+tensor2fourindex(t,ff);
+cout <<"Error = "<<(f-ff).norm()<<endl;
+}
+
+
+
 
 }//main

tensor.h

@@ -266,22 +266,26 @@ public:
 
 
 //due to template nesting and specialization limitation, this cannot be class member - the S parameter has to be the outer one
+//note also that the class Tnesor does not support symmetry between groups of indices (only symmetry inside each group), the full fourindex symmetry cannot thus be reflected in some cases
+//
 template<fourindexsymtype S, typename T, typename I> 
 Tensor<T> fourindex2tensor(const fourindex_dense<S,T,I> &f);
 template<fourindexsymtype S, typename T, typename I>          
 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
 template<typename T, typename I>
 Tensor<T>  fourindex2tensor(const fourindex_dense<nosymmetry,T,I> &f)
 {
 NRVec<indexgroup> shape(4);
+int n=f.nbas();
 for(int i=0; i<4; ++i)
 	{
 	shape[i].number=1;
-	shape[i].symmetry-0;
+	shape[i].symmetry=0;
-	shape[i].offset=0;
+	shape[i].offset=1;
-	shape[i].range=f.nbas();
+	shape[i].range=n;
 	}
 NRVec<T> data(f);
 return Tensor<T>(shape,data);
@@ -303,6 +307,114 @@ f=fourindex_dense<nosymmetry,T,I>(range,NRMat<T>(t.data,range*range,range*range)
 }
 
 
+template<typename T, typename I>
+Tensor<T>  fourindex2tensor(const fourindex_dense<twoelectronrealmullikanAB,T,I> &f)
+{
+NRVec<indexgroup> shape(2);
+int n=f.nbas();
+for(int i=0; i<2; ++i)
+	{
+	shape[i].number=2;
+	shape[i].symmetry=1;
+	shape[i].offset=1;
+	shape[i].range=n;
+	}
+NRVec<T> data(f);
+return Tensor<T>(shape,data);
+}
+
+template<typename T, typename I>
+void tensor2fourindex(const Tensor<T> &t, fourindex_dense<twoelectronrealmullikanAB,T,I> &f)
+{
+if(t.rank()!=4) laerror("wrong rank in tensor2fourindex");
+if(t.shape.size()!=2) laerror("wrong symmetry groups in tensor2fourindex");
+int range=t.shape[0].range;
+int offset=t.shape[0].offset;
+for(int i=0; i<t.shape.size(); ++i)
+	{
+	if(range!=t.shape[i].range) laerror("range mismatch in tensor2fourindex");
+	if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
+	if(1!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
+	}
+f=fourindex_dense<twoelectronrealmullikanAB,T,I>(NRMat<T>(t.data,range*(range+1)/2,range*(range+1)/2));
+}
+
+
+template<typename T, typename I>
+Tensor<T>  fourindex2tensor(const fourindex_dense<twoelectronrealmullikan,T,I> &f)
+{
+NRVec<indexgroup> shape(2);
+int n=f.nbas();
+for(int i=0; i<2; ++i)
+	{
+	shape[i].number=2;
+	shape[i].symmetry=1;
+	shape[i].offset=1;
+	shape[i].range=n;
+	}
+NRMat<T> mat(f);
+NRVec<T> data(mat); //expand symmetric to full matrix for index group symmetry
+return Tensor<T>(shape,data);
+}
+
+template<typename T, typename I>
+void tensor2fourindex(const Tensor<T> &t, fourindex_dense<twoelectronrealmullikan,T,I> &f)
+{
+if(t.rank()!=4) laerror("wrong rank in tensor2fourindex");
+if(t.shape.size()!=2) laerror("wrong symmetry groups in tensor2fourindex");
+int range=t.shape[0].range;
+int offset=t.shape[0].offset;
+for(int i=0; i<t.shape.size(); ++i)
+	{
+	if(range!=t.shape[i].range) laerror("range mismatch in tensor2fourindex");
+	if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
+	if(1!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
+	}
+NRMat<T> mat(t.data,range*(range+1)/2,range*(range+1)/2);
+f=fourindex_dense<twoelectronrealmullikan,T,I>(NRSMat<T>(mat)); //symmetrize mat 
+}
+
+
+template<typename T, typename I>
+Tensor<T>  fourindex2tensor(const fourindex_dense<T2IjAb_aces,T,I> &f)
+{
+NRVec<indexgroup> shape(4);
+for(int i=0; i<4; ++i)
+	{
+	shape[i].number=1;
+	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;
+NRVec<T> data(f);
+return Tensor<T>(shape,data);
+}
+
+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 offset=t.shape[0].offset;
+for(int i=0; i<t.shape.size(); ++i)
+	{
+	if(offset!=t.shape[i].offset) laerror("offset mismatch in tensor2fourindex");
+	if(0!=t.shape[i].symmetry) laerror("symmetry mismatch in tensor2fourindex");
+	}
+NRMat<T> mat(t.data,noca*nocb,nvra*nvrb);
+f=fourindex_dense<T2IjAb_aces,T,I>(noca,nocb,nvra,nvrb,mat);
+}
+
+
+
+