Crawdad 双电子积分维数转换问题

最爱喵星人

我想把存为一维数组的双电子积分转成四维矩阵。（一维存储双电子积分：http://sirius.chem.vt.edu/~crawd ... 3/h2o_sto3g/eri.dat）

根据下面链接里step #3提到的方法
http://sirius.chem.vt.edu/wiki/d ... rogramming:project3

我用下面的代码输出双电子积分的index,但是却与http://sirius.chem.vt.edu/~crawd ... 3/h2o_sto3g/eri.dat 差很多。请问哪里有问题？谢谢。

#include<iostream>
using namespace std;
int main() {

   int i,j,k,l,ij,kl,ijkl;
   for(i=0; i<7; i++){
     for(j=0; j<=i; j++){
       ij=i*(i+1)/2+j;
       for(k=0; k<7; k++){
          for(l=0; l<=k;l++){
             kl=k*(k+1)/2+l;
             if(ij>=kl)
             {
             ijkl=ij*(ij+1)/2+kl;
             cout<<ijkl<<" "<<i+1<<" "<<j+1<<" "<<k+1<<" "<<l+1;
             }
             cout<<endl;
     }
   }
  }
}

return 0;
}

Warm_Cloud

双电子积分实际上是一个shell一个shell算，然后每个shell包含不同的bas

Warm_Cloud

    ii = 0
    do ipsh = 1,ShlPar
        ish = Pos_par(1,ipsh)
        jsh = Pos_par(2,ipsh)
        kbas(1) = PosShl(ish)
        kbas(2) = PosShl(jsh)
        mbas(1) = shl_num(ish)
        mbas(2) = shl_num(jsh)
      
        do jpsh = 1,ShlPar
            ksh = Pos_par(1,jpsh)
            lsh = Pos_par(2,jpsh)
            kbas(3) = PosShl(ksh)
            kbas(4) = PosShl(lsh)
            mbas(3) = shl_num(ksh)
            mbas(4) = shl_num(lsh)
               
            if(IJud1(ipsh) < IJud2(jpsh) ) cycle      
            if( TE(ipsh)*TE(jpsh) < intthre ) cycle

            do i = 1,mbas(1)
               i1 = kbas(1)+i-1
               do j = 1,mbas(2)
                  j1 = kbas(2)+j-1
                  do k = 1,mbas(3)
                     k1 = kbas(3)+k-1
                     do l = 1,mbas(4)
                        l1 = kbas(4)+l-1
                           
                        if( i1 < j1 ) cycle
                        if( k1 < l1 ) cycle
                        if( i1*(i1-1)+2*j1 < k1*(k1-1)+2*l1 ) cycle
                        ii = ii + 1
                        Temp1 = Two_Ele(ii)
                        
                        g_Mo(EI(i1,j1),EI(k1,l1)) = Temp1
                        g_Mo(EI(k1,l1),EI(i1,j1)) = Temp1

                     end do
                  end do
               end do
            end do
               
       end do
    end do