From 04b4ef9f72515eb0ec1cf067292e56f5cad6dda3 Mon Sep 17 00:00:00 2001 From: jiri Date: Sat, 9 Sep 2006 19:51:04 +0000 Subject: [PATCH] *** empty log message *** --- fourindex.h | 96 ++++++++++++++++++++++++++++++++++++++++++++--------- 1 file changed, 81 insertions(+), 15 deletions(-) diff --git a/fourindex.h b/fourindex.h index ff87971..ec6dbad 100644 --- a/fourindex.h +++ b/fourindex.h @@ -78,9 +78,9 @@ struct matel4stored //later add symmetry of complex integrals -typedef enum {undefined_symmetry=-1,nosymmetry=0, twoelectronrealmullikan=1, twoelectronrealdirac=2, T2ijab_real=3, antisymtwoelectronrealmullikan=4,antisymtwoelectronrealdirac=5} fourindexsymtype; //only permutation-nonequivalent elements are stored +typedef enum {undefined_symmetry=-1,nosymmetry=0, twoelectronrealmullikan=1, twoelectronrealdirac=2, T2ijab_aces=3, antisymtwoelectronrealmullikan=4,antisymtwoelectronrealdirac=5, T2IjAb_aces=6} fourindexsymtype; //only permutation-nonequivalent elements are stored // these should actually be static private members of the fourindex class, but leads to an ICE on gcc3.2 -static const int fourindex_n_symmetrytypes=6; +static const int fourindex_n_symmetrytypes=7; static const int fourindex_permnumbers[fourindex_n_symmetrytypes]={1,8,8,4,16,16}; static const int fourindex_permutations[fourindex_n_symmetrytypes][16][5]= { @@ -90,6 +90,7 @@ static const int fourindex_permutations[fourindex_n_symmetrytypes][16][5]= {{0,1,2,3,1},{1,0,2,3,-1},{0,1,3,2,-1},{1,0,3,2,1}}, {{0,1,2,3,1}, {1,0,2,3,1}, {0,1,3,2,1}, {1,0,3,2,1}, {2,3,0,1,1}, {3,2,0,1,1}, {2,3,1,0,1}, {3,2,1,0,1}, {0,3,2,1,-1}, {1,3,2,0,-1}, {0,2,3,1,-1}, {1,2,3,1,-1}, {2,1,0,3,-1}, {3,1,0,2,-1}, {2,0,1,3,-1}, {3,0,1,2,-1}}, {{0,1,2,3,1},{2,1,0,3,1},{0,3,2,1,1},{2,3,0,1,1},{1,0,3,2,1},{3,0,1,2,1},{1,2,3,0,1},{3,2,1,0,1}, {0,1,3,2,-1},{2,1,3,0,-1},{0,3,1,2,-1},{2,3,1,0,-1},{1,0,2,3,-1},{3,0,2,1,-1},{1,2,0,3,-1},{3,2,0,1,-1}}, + {{0,1,2,3,1}}, //like nosymmetry but different index ranges }; @@ -112,7 +113,8 @@ switch(symmetry) if(index.indiv.i==index.indiv.j && index.indiv.k==index.indiv.l || index.indiv.k==index.indiv.j && index.indiv.i==index.indiv.l) elem*=.5; break; - case T2ijab_real: break; //result will automatically vanish due to generated antisymmetry; i!=a from principle + case T2ijab_aces: break; //result will automatically vanish due to generated antisymmetry; i!=a from principle + case T2IjAb_aces: break; //no actual symmetry case nosymmetry: break; default: laerror("illegal symmetry"); } @@ -209,6 +211,10 @@ public: inline void add(const I (&index)[4], const T elem) {matel4 *ltmp= new matel4; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &index, sizeof(union packed_index)); list->elem=elem;} + inline void add(const matel4 &rhs) + {matel4 *ltmp= new matel4; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &rhs.index, sizeof(union packed_index)); list->elem=rhs.elem;} + inline void add(const matel4stored &rhs) + {matel4 *ltmp= new matel4; ltmp->next=list; list=ltmp; memcpy(&list->index.packed, &rhs.index, sizeof(union packed_index)); list->elem=rhs.elem;} unsigned long put(int fd,bool withattr=true) const; unsigned long get(int fd,bool withattr=true); }; @@ -585,6 +591,7 @@ istream& operator>>(istream &s, fourindex &x) return s; } + /////////////////////densely stored fourindex/////////////////////////////////// //not all symmetry cases implemented yet, but a general template declaration used //we use a general template forward declaration, but then it has to be done differently for (almost) each case @@ -620,9 +627,9 @@ fourindex_dense::fourindex_dense::iterator p; #ifdef DEBUG -unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j); -unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l); -if (I<0 || I>=(unsigned long)NRSMat::nn || J<0 || J>=(unsigned long)NRSMat::nn) laerror("fourindex_dense index out of range in constructor"); +unsigned int IJ = SMat_index_1(p->index.indiv.i,p->index.indiv.j); +unsigned int KL = SMat_index_1(p->index.indiv.k,p->index.indiv.l); +if (IJ<0 || IJ>=(unsigned int)NRSMat::nn || KL<0 || KL>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range in constructor"); #endif for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j,p->index.indiv.k,p->index.indiv.l) = p->elem; } @@ -633,9 +640,9 @@ fourindex_dense::fourindex_dense::iterator p; #ifdef DEBUG -unsigned long I = SMat_index_1(p->index.indiv.i,p->index.indiv.j); -unsigned long J = SMat_index_1(p->index.indiv.k,p->index.indiv.l); -if (I<0 || I>=(unsigned long)NRSMat::nn || J<0 || J>=(unsigned long)NRSMat::nn) laerror("fourindex_dense index out of range in constructor"); +unsigned int IJ = SMat_index_1(p->index.indiv.i,p->index.indiv.j); +unsigned int KL = SMat_index_1(p->index.indiv.k,p->index.indiv.l); +if (IJ<0 || IJ>=(unsigned int)NRSMat::nn || KL<0 || KL>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range in constructor"); #endif for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j ,p->index.indiv.k,p->index.indiv.l) = p->elem; } @@ -645,12 +652,12 @@ for(p=rhs.begin(); p!= rhs.end(); ++p) (*this)(p->index.indiv.i,p->index.indiv.j template T& fourindex_dense::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) { -unsigned long I = SMat_index_1(i,j); -unsigned long J = SMat_index_1(k,l); +unsigned int I = SMat_index_1(i,j); +unsigned int J = SMat_index_1(k,l); //I,J act as indices of a NRSmat #ifdef DEBUG if (*NRSMat::count != 1) laerror("lval (i,j,k,l) with count > 1 in fourindex_dense"); - if (I<0 || I>=(unsigned long)NRSMat::nn || J<0 || J>=(unsigned long)NRSMat::nn) laerror("fourindex_dense index out of range"); + if (I<0 || I>=(unsigned int)NRSMat::nn || J<0 || J>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range"); if (!NRSMat::v) laerror("access to unallocated fourindex_dense"); #endif return NRSMat::v[SMat_index(I,J)]; @@ -660,17 +667,76 @@ return NRSMat::v[SMat_index(I,J)]; template const T& fourindex_dense::operator() (unsigned int i, unsigned int j, unsigned int k, unsigned int l) const { -unsigned long I = SMat_index_1(i,j); -unsigned long J = SMat_index_1(k,l); +unsigned int I = SMat_index_1(i,j); +unsigned int J = SMat_index_1(k,l); //I,J act as indices of a NRSmat #ifdef DEBUG - if (I<0 || I>=(unsigned long)NRSMat::nn || J<0 || J>=(unsigned long)NRSMat::nn) laerror("fourindex_dense index out of range"); + if (I<0 || I>=(unsigned int)NRSMat::nn || J<0 || J>=(unsigned int)NRSMat::nn) laerror("fourindex_dense index out of range"); if (!NRSMat::v) laerror("access to unallocated fourindex_dense"); #endif return NRSMat::v[SMat_index(I,J)]; } +//access to spin-blocks of T2 amplitudes in aces storage order +//both occupied and virtual indices start from 1 + +template +class fourindex_dense : public NRMat { +private: + unsigned int noca,nocb,nvra,nvrb; +public: + fourindex_dense(): NRMat() {noca=nocb=nvra=nvrb=0;}; + explicit fourindex_dense(const int nocca, const int noccb, const int nvrta, const int nvrtb): NRMat(nocca*noccb,nvrta*nvrtb) {noca=nocca; nocb=noccb; nvra=nvrta; nvrb=nvrtb;}; + +//here i,a are alpha j,b beta + inline T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) +{ +#ifdef DEBUG +if(i<1||i>noca ||j<1||j>nocb|| a<1||a>nvra||b<1||b>nvrb) laerror("T2IjAb_aces fourindex out of range") +if (!NRMat::v) laerror("access to unallocated fourindex_dense"); +#endif +return (*this).NRMat::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1); +} + inline const T& operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) const +{ +#ifdef DEBUG +if(i<1||i>noca ||j<1||j>nocb|| a<1||a>nvra||b<1||b>nvrb) laerror("T2IjAb_aces fourindex out of range") +if (!NRMat::v) laerror("access to unallocated fourindex_dense"); +#endif +return (*this).NRMat::operator() ((j-1)*noca+i-1,(b-1)*nvra+a-1); +} +}; + + +template +class fourindex_dense : public NRMat { +private: + unsigned int nocc,nvrt,ntri; +public: + fourindex_dense(): NRMat() {nocc=nvrt=ntri=0;}; + explicit fourindex_dense(const int noc, const int nvr): NRMat(noc*(noc-1)/2,nvr*(nvr-1)/2) {nocc=noc; nvrt=nvr; ntri=nvr*(nvr-1)/2;}; + +//we cannot return reference due to the possible sign change +//stored values are for i>j a>b + inline T operator() (unsigned int i, unsigned int j, unsigned int a, unsigned int b) +{ +#ifdef DEBUG +if(i<1||i>nocc ||j<1||j>nocc|| a<1||a>nvrt||b<1||b>nvrt) laerror("T2ijab_aces fourindex out of range") +if (!NRMat::v) laerror("access to unallocated fourindex_dense"); +#endif +int minus=0; +if(i::operator() ((i-2)*(i-1)/2+j-1, (a-2)*(a-1)/2+b-1); +return (minus&1)? -val:val; +} + +}; + + + + #endif /*_fourindex_included*/