i依米花 2020-10-14 18:40
浏览 40

最近在dev c++写了一个用rk4模拟自旋系统运动的代码，出现了一个奇怪的错误，为什么代码有时候可以运行，有时候会卡住？

这个代码有时候会运行成功，但是重复运行的时候偶尔又会出错。怀疑是内存溢出，但是溢出这个问题的出现会是随机的吗？

#include
#include
#include
#include
#include
#include
#include
//system parameters
#define L 8/* the size of the system,total number of spin=L^3*4 */
#define V 2048
#define NB 6
#define PI 3.1415926535
#define SAMPLE_NUMBER 1024
#define step 8192//自旋状态更新次数

#define MAX_LEN 256

//#define T /*the temperature */

static double T=0.01; //温度
static double B=0.0; //外场
static double t=0.0; //时间

static double dT=0.01;
static double dB=0.1;

static double steplength=0.025;//spin运动的时间间隔

static double state[3][V]={0};

static int nn[NB][V]={0}; //计算每个spin的最近邻自旋

int number (int r1,int r2,int r3,int r4)
{
return 4*L*L*r1+4*L*r2+4*r3+r4;
}

void near (int r1,int r2,int r3,int r4,int neighbor[NB])
{

if (r4==0)
{
    neighbor[0]=number(r1,r2,r3,1);
    neighbor[1]=number(r1,r2,r3,2);
    neighbor[2]=number(r1,r2,r3,3);
    neighbor[3]=number((r1-1+L)%L,r2,r3,1);
    neighbor[4]=number(r1,(r2-1+L)%L,r3,2);
    neighbor[5]=number(r1,r2,(r3-1+L)%L,3); 
}

else if(r4==1)  
{
    neighbor[0]=number(r1,r2,r3,0);
    neighbor[1]=number(r1,r2,r3,2);
    neighbor[2]=number(r1,r2,r3,3);
    neighbor[3]=number((r1+1)%L,r2,r3,0);
    neighbor[4]=number((r1+1)%L,(r2-1+L)%L,r3,2);
    neighbor[5]=number((r1+1)%L,r2,(r3-1+L)%L,3);   
}

else if(r4==2)
{
    neighbor[0]=number(r1,r2,r3,0);
    neighbor[1]=number(r1,r2,r3,1);
    neighbor[2]=number(r1,r2,r3,3);
    neighbor[3]=number(r1,(r2+1)%L,r3,0);
    neighbor[4]=number((r1-1+L)%L,(r2+1)%L,r3,1);
    neighbor[5]=number(r1,(r2+1)%L,(r3-1+L)%L,3);   

}

else if(r4==3)
{                                                                                                                                                                                                                
    neighbor[0]=number(r1,r2,r3,0);
    neighbor[1]=number(r1,r2,r3,1);
    neighbor[2]=number(r1,r2,r3,2);
    neighbor[3]=number(r1,r2,(r3+1)%L,0);
    neighbor[4]=number((r1-1+L)%L,r2,(r3+1)%L,1);
    neighbor[5]=number(r1,(r2-1+L)%L,(r3+1)%L,2);           
}

}

//calculate spin_structure_factor
void table1 (int cn[4][V])
{

int r1=0,r2=0,r3=0,r4=0,i=0;
for(r1=0;r1<L;r1++)
{
    for(r2=0;r2<L;r2++)
    {
        for(r3=0;r3<L;r3++)
        {
            for(r4=0;r4<4;r4++)
            {
                i=number(r1,r2,r3,r4);
                cn[3][i]=r4;cn[2][i]=r3;cn[1][i]=r2;cn[0][i]=r1;


            }
        }
    }
}

}

//for a special wave vector,get the phase of every spin
//wave vector=2PI*[hhk],h=integer/L;

void phase(double wvector[3],double real[V],double image[V],int cn[4][V])
{
int i=0;

double angle=0,wvector1[3]={0};

wvector1[0]=1.0*(wvector[0]+wvector[2])/2.0;
wvector1[1]=1.0*(wvector[0]+wvector[1])/2.0;
wvector1[2]=1.0*(wvector[1]+wvector[2])/2.0;

//计算K点以b1=(1,-1,1),b2=(1,1,-1),b3=(-1,1,1)为基矢的坐标 

for(i=0;i<V;i++)
{
    if(cn[3][i]==0)
    {
        angle=2*PI*(wvector1[0]*(cn[0][i]-0.125)+wvector1[1]*(cn[1][i]-0.125)+wvector1[2]*(cn[2][i]-0.125));
        real[i]=cos(angle);
        image[i]=sin(angle);
    }
    else if(cn[3][i]==1)
    {
        angle=2*PI*(wvector1[0]*(cn[0][i]+0.375)+wvector1[1]*(cn[1][i]-0.125)+wvector1[2]*(cn[2][i]-0.125));
        real[i]=cos(angle);
        image[i]=sin(angle);    
    }
    else if(cn[3][i]==2)
    {
        angle=2*PI*(wvector1[0]*(cn[0][i]-0.125)+wvector1[1]*(cn[1][i]+0.375)+wvector1[2]*(cn[2][i]-0.125));
        real[i]=cos(angle);
        image[i]=sin(angle);    

    }
    else if(cn[3][i]==3)    
    {
        angle=2*PI*(wvector1[0]*(cn[0][i]-0.125)+wvector1[1]*(cn[1][i]-0.125)+wvector1[2]*(cn[2][i]+0.375));
        real[i]=cos(angle);
        image[i]=sin(angle);    
    }       
}

}

void cal_sq(double s[3][V],double real[V],double image[V],double sq[3][2])
{
int i=0,j=0;
for(i=0;i<3;i++)
{
for(j=0;j<2;j++)
{
sq[i][j]=0;

}

}
//printf("state[0][0]=%lf\n",s[0][0]);printf("state[1][0]=%lf\n",s[1][0]);

for(j=0;j<3;j++)/*指标j是对结构因子的三个分量x,y,z进行计算*/ 
{   
    for(i=0;i<V;i++)
    {
        sq[j][0]=sq[j][0]+1.0*s[j][i]*real[i];
        sq[j][1]=sq[j][1]+1.0*s[j][i]*image[i];
        //printf("%d    %lf   %lf   %lf   %lf   %lf   %lf\n",i,sqr[0],sqi[0],sqr[1],sqi[1],sqr[2],sqi[2]);
    }
}

}

//计算结构因子的实部和虚部
void cal_dssf(double sq_0[3][2], double sq_1[3][2],double dssf[2])
{
int i=0;

for(i=0;i<3;i++)
{
    dssf[0]=dssf[0]+(sq_0[i][0]*sq_1[i][0]-sq_0[i][1]*sq_1[i][1]);
    dssf[1]=dssf[1]+(sq_0[i][0]*sq_1[i][1]+sq_0[i][1]*sq_1[i][0]);
} 
dssf[0]=1.0*dssf[0]/(L*L*L);
dssf[1]=1.0*dssf[1]/(L*L*L);

}

void Complex_conjugation(double sq_0[3][2],double sq_1[3][2]) //取复共轭，三维的复矢量
{
int i=0;
for(i=0;i<3;i++)
{
sq_0[i][0]=sq_1[i][0];
sq_0[i][1]=-sq_1[i][1];
}

}

//nn[N][B]保存了所有自旋的最近邻自旋
void table (int nn[NB][V])
{
int r1=0,r2=0,r3=0,r4=0,index=0,i=0;
int neighbor[NB]={0};

for(r1=0;r1<L;r1++)
{
    for(r2=0;r2<L;r2++)
    {
        for(r3=0;r3<L;r3++)
        {
            for(r4=0;r4<4;r4++)
            {
               index=number(r1,r2,r3,r4);
               near(r1,r2,r3,r4,neighbor);

               for(i=0;i<NB;i++)
               {
                  nn[i][index]=neighbor[i];
               }

            }

        }
    }

}

}

//计算一个spin在有效的磁场作用下的运动
double vector_length (double vector[3])
{
return sqrt(vector[0]*vector[0]+vector[1]*vector[1]+vector[2]*vector[2]);
}

void spin_move(double s[3][V],double B)
{

 double k1[3][V]={0},k2[3][V]={0},k3[3][V]={0},k4[3][V]={0};

 int i=0,j=0,n[6]={0};

 //分别计算RK4的四个系数 

      for(i=0;i<V;i++)
   {

      for(j=0;j<6;j++)   //对于每一个自旋，首先找到它的最近邻自旋 
      {
        n[j]=nn[j][i];                  
      }

     //fprintf(fpn,"%d      %d       %d       %d      %d      %d        %d\n",i,n[0],n[1],n[2],n[3],n[4],n[5]); 
     k1[0][i]=-(s[1][i]*(s[2][n[0]]+s[2][n[1]]+s[2][n[2]]+s[2][n[3]]+s[2][n[4]]+s[2][n[5]]+B)-s[2][i]*(s[1][n[0]]+s[1][n[1]]+s[1][n[2]]+s[1][n[3]]+s[1][n[4]]+s[1][n[5]]));
     k1[1][i]=-(s[2][i]*(s[0][n[0]]+s[0][n[1]]+s[0][n[2]]+s[0][n[3]]+s[0][n[4]]+s[0][n[5]])-s[0][i]*(s[2][n[0]]+s[2][n[1]]+s[2][n[2]]+s[2][n[3]]+s[2][n[4]]+s[2][n[5]]+B));
     k1[2][i]=-(s[0][i]*(s[1][n[0]]+s[1][n[1]]+s[1][n[2]]+s[1][n[3]]+s[1][n[4]]+s[1][n[5]])-s[1][i]*(s[0][n[0]]+s[0][n[1]]+s[0][n[2]]+s[0][n[3]]+s[0][n[4]]+s[0][n[5]]));

        //fprintf(fpk1,"%d      %12.8lf       %12.8lf      %12.8lf      \n",i,k1[0][i],k1[1][i],k1[2][i]);  
   }   

      for(i=0;i<V;i++)
   {

      for(j=0;j<6;j++)   //对于每一个自旋，首先找到它的最近邻自旋 
      {
        n[j]=nn[j][i];                  
      }

     //fprintf(fpn,"%d      %d       %d       %d      %d      %d        %d\n",i,n[0],n[1],n[2],n[3],n[4],n[5]); 
     k2[0][i]=-((s[1][i]+0.5*steplength*k1[1][i])*((s[2][n[0]]+0.5*steplength*k1[2][n[0]])+(s[2][n[1]]+0.5*steplength*k1[2][n[1]])+(s[2][n[2]]+0.5*steplength*k1[2][n[2]])+(s[2][n[3]]+0.5*steplength*k1[2][n[3]])+(s[2][n[4]]+0.5*steplength*k1[2][n[4]])+(s[2][n[5]]+0.5*steplength*k1[2][n[5]])+B)-(s[2][i]+0.5*steplength*k1[2][i])*((s[1][n[0]]+0.5*steplength*k1[1][n[0]])+(s[1][n[1]]+0.5*steplength*k1[1][n[1]])+(s[1][n[2]]+0.5*steplength*k1[1][n[2]])+(s[1][n[3]]+0.5*steplength*k1[1][n[3]])+(s[1][n[4]]+0.5*steplength*k1[1][n[4]])+(s[1][n[5]]+0.5*steplength*k1[1][n[5]])));
     k2[1][i]=-((s[2][i]+0.5*steplength*k1[2][i])*((s[0][n[0]]+0.5*steplength*k1[0][n[0]])+(s[0][n[1]]+0.5*steplength*k1[0][n[1]])+(s[0][n[2]]+0.5*steplength*k1[0][n[2]])+(s[0][n[3]]+0.5*steplength*k1[0][n[3]])+(s[0][n[4]]+0.5*steplength*k1[0][n[4]])+(s[0][n[5]]+0.5*steplength*k1[0][n[5]]))-(s[0][i]+0.5*steplength*k1[0][i])*((s[2][n[0]]+0.5*steplength*k1[2][n[0]])+(s[2][n[1]]+0.5*steplength*k1[2][n[1]])+(s[2][n[2]]+0.5*steplength*k1[2][n[2]])+(s[2][n[3]]+0.5*steplength*k1[2][n[3]])+(s[2][n[4]]+0.5*steplength*k1[2][n[4]])+(s[2][n[5]]+0.5*steplength*k1[2][n[5]])+B));
     k2[2][i]=-((s[0][i]+0.5*steplength*k1[0][i])*((s[1][n[0]]+0.5*steplength*k1[1][n[0]])+(s[1][n[1]]+0.5*steplength*k1[1][n[1]])+(s[1][n[2]]+0.5*steplength*k1[1][n[2]])+(s[1][n[3]]+0.5*steplength*k1[1][n[3]])+(s[1][n[4]]+0.5*steplength*k1[1][n[4]])+(s[1][n[5]]+0.5*steplength*k1[1][n[5]]))-(s[1][i]+0.5*steplength*k1[1][i])*((s[0][n[0]]+0.5*steplength*k1[0][n[0]])+(s[0][n[1]]+0.5*steplength*k1[0][n[1]])+(s[0][n[2]]+0.5*steplength*k1[0][n[2]])+(s[0][n[3]]+0.5*steplength*k1[0][n[3]])+(s[0][n[4]]+0.5*steplength*k1[0][n[4]])+(s[0][n[5]]+0.5*steplength*k1[0][n[5]])));



        //fprintf(fpk2,"%d      %12.8lf       %12.8lf      %12.8lf      \n",i,k2[0][i],k2[1][i],k2[2][i]);  
   }       

      for(i=0;i<V;i++)
   {

      for(j=0;j<6;j++)   //对于每一个自旋，首先找到它的最近邻自旋 
      {
        n[j]=nn[j][i];                  
      }

     //fprintf(fpn,"%d      %d       %d       %d      %d      %d        %d\n",i,n[0],n[1],n[2],n[3],n[4],n[5]); 
     k3[0][i]=-((s[1][i]+0.5*steplength*k2[1][i])*((s[2][n[0]]+0.5*steplength*k2[2][n[0]])+(s[2][n[1]]+0.5*steplength*k2[2][n[1]])+(s[2][n[2]]+0.5*steplength*k2[2][n[2]])+(s[2][n[3]]+0.5*steplength*k2[2][n[3]])+(s[2][n[4]]+0.5*steplength*k2[2][n[4]])+(s[2][n[5]]+0.5*steplength*k2[2][n[5]])+B)-(s[2][i]+0.5*steplength*k2[2][i])*((s[1][n[0]]+0.5*steplength*k2[1][n[0]])+(s[1][n[1]]+0.5*steplength*k2[1][n[1]])+(s[1][n[2]]+0.5*steplength*k2[1][n[2]])+(s[1][n[3]]+0.5*steplength*k2[1][n[3]])+(s[1][n[4]]+0.5*steplength*k2[1][n[4]])+(s[1][n[5]]+0.5*steplength*k2[1][n[5]])));
     k3[1][i]=-((s[2][i]+0.5*steplength*k2[2][i])*((s[0][n[0]]+0.5*steplength*k2[0][n[0]])+(s[0][n[1]]+0.5*steplength*k2[0][n[1]])+(s[0][n[2]]+0.5*steplength*k2[0][n[2]])+(s[0][n[3]]+0.5*steplength*k2[0][n[3]])+(s[0][n[4]]+0.5*steplength*k2[0][n[4]])+(s[0][n[5]]+0.5*steplength*k2[0][n[5]]))-(s[0][i]+0.5*steplength*k2[0][i])*((s[2][n[0]]+0.5*steplength*k2[2][n[0]])+(s[2][n[1]]+0.5*steplength*k2[2][n[1]])+(s[2][n[2]]+0.5*steplength*k2[2][n[2]])+(s[2][n[3]]+0.5*steplength*k2[2][n[3]])+(s[2][n[4]]+0.5*steplength*k2[2][n[4]])+(s[2][n[5]]+0.5*steplength*k2[2][n[5]])+B));
     k3[2][i]=-((s[0][i]+0.5*steplength*k2[0][i])*((s[1][n[0]]+0.5*steplength*k2[1][n[0]])+(s[1][n[1]]+0.5*steplength*k2[1][n[1]])+(s[1][n[2]]+0.5*steplength*k2[1][n[2]])+(s[1][n[3]]+0.5*steplength*k2[1][n[3]])+(s[1][n[4]]+0.5*steplength*k2[1][n[4]])+(s[1][n[5]]+0.5*steplength*k2[1][n[5]]))-(s[1][i]+0.5*steplength*k2[1][i])*((s[0][n[0]]+0.5*steplength*k2[0][n[0]])+(s[0][n[1]]+0.5*steplength*k2[0][n[1]])+(s[0][n[2]]+0.5*steplength*k2[0][n[2]])+(s[0][n[3]]+0.5*steplength*k2[0][n[3]])+(s[0][n[4]]+0.5*steplength*k2[0][n[4]])+(s[0][n[5]]+0.5*steplength*k2[0][n[5]])));



        //fprintf(fpk3,"%d      %12.8lf       %12.8lf      %12.8lf      \n",i,k3[0][i],k3[1][i],k3[2][i]);  
   }       

      for(i=0;i<V;i++)
   {

      for(j=0;j<6;j++)   //对于每一个自旋，首先找到它的最近邻自旋 
      {
        n[j]=nn[j][i];                  
      }

     //fprintf(fpn,"%d      %d       %d       %d      %d      %d        %d\n",i,n[0],n[1],n[2],n[3],n[4],n[5]); 
     k4[0][i]=-((s[1][i]+steplength*k3[1][i])*((s[2][n[0]]+steplength*k3[2][n[0]])+(s[2][n[1]]+steplength*k3[2][n[1]])+(s[2][n[2]]+steplength*k3[2][n[2]])+(s[2][n[3]]+steplength*k3[2][n[3]])+(s[2][n[4]]+steplength*k3[2][n[4]])+(s[2][n[5]]+steplength*k3[2][n[5]])+B)-(s[2][i]+steplength*k3[2][i])*((s[1][n[0]]+steplength*k3[1][n[0]])+(s[1][n[1]]+steplength*k3[1][n[1]])+(s[1][n[2]]+steplength*k3[1][n[2]])+(s[1][n[3]]+steplength*k3[1][n[3]])+(s[1][n[4]]+steplength*k3[1][n[4]])+(s[1][n[5]]+steplength*k3[1][n[5]])));
     k4[1][i]=-((s[2][i]+steplength*k3[2][i])*((s[0][n[0]]+steplength*k3[0][n[0]])+(s[0][n[1]]+steplength*k3[0][n[1]])+(s[0][n[2]]+steplength*k3[0][n[2]])+(s[0][n[3]]+steplength*k3[0][n[3]])+(s[0][n[4]]+steplength*k3[0][n[4]])+(s[0][n[5]]+steplength*k3[0][n[5]]))-(s[0][i]+steplength*k3[0][i])*((s[2][n[0]]+steplength*k3[2][n[0]])+(s[2][n[1]]+steplength*k3[2][n[1]])+(s[2][n[2]]+steplength*k3[2][n[2]])+(s[2][n[3]]+steplength*k3[2][n[3]])+(s[2][n[4]]+steplength*k3[2][n[4]])+(s[2][n[5]]+steplength*k3[2][n[5]])+B));
     k4[2][i]=-((s[0][i]+steplength*k3[0][i])*((s[1][n[0]]+steplength*k3[1][n[0]])+(s[1][n[1]]+steplength*k3[1][n[1]])+(s[1][n[2]]+steplength*k3[1][n[2]])+(s[1][n[3]]+steplength*k3[1][n[3]])+(s[1][n[4]]+steplength*k3[1][n[4]])+(s[1][n[5]]+steplength*k3[1][n[5]]))-(s[1][i]+steplength*k3[1][i])*((s[0][n[0]]+steplength*k3[0][n[0]])+(s[0][n[1]]+steplength*k3[0][n[1]])+(s[0][n[2]]+steplength*k3[0][n[2]])+(s[0][n[3]]+steplength*k3[0][n[3]])+(s[0][n[4]]+steplength*k3[0][n[4]])+(s[0][n[5]]+steplength*k3[0][n[5]])));
        //fprintf(fpk4,"%d      %12.8lf       %12.8lf      %12.8lf      \n",i,k4[0][i],k4[1][i],k4[2][i]);  
   }       

    for(i=0;i<V;i++)
   {
       for(j=0;j<3;j++)
      {
       s[j][i]=s[j][i]+1.0*steplength*(k1[j][i]+2.0*k2[j][i]+2.0*k3[j][i]+k4[j][i])/6.0;
       }

    }

}

void set_zero(double sq_0[3][2])
{
int i=0,j=0;
for(i=0;i<3;i++)
{
for(j=0;j<2;j++)
{
sq_0[i][j]=0;
}
}

}

double tote (double s[3][V],int nb[NB][V],double ext_field)
{

int i=0,j=0,n=0;
double sum=0;

for(i=0;i<V;i++)
{
    for(j=0;j<NB;j++)
    {
        n=nb[j][i];
        //printf("%5d   ",n);
        sum=sum+s[0][i]*s[0][n]+s[1][i]*s[1][n]+s[2][i]*s[2][n];
    }
   // printf("\n");
}

sum=sum*0.5;

for(i=0;i<V;i++)
{
    sum=sum-s[2][i]*ext_field;
}

return sum/(1.0*V);

//return 0.5*sum/(1.0*V);

}

void mom (double s[3][V],double a[3])
{
int i=0,j=0;

for(j=0;j<3;j++)
{
    for(i=0;i<V;i++)
    {
        a[j]=a[j]+1.0*s[j][i];
    }

}

}

double dot_product(double a[3] ,double b[3])
{
return a[0]*b[0]+a[1]*b[1]+a[2]*b[2];

}

int main()
{
double b[3][V]={0},state[3][V]={0},t=0;
//double s0[3]={0},s1[3]={0};//s0保存了初始
double sq_0[3][2]={0},sq_1[3][2]={0};//保存了S(0,-q)和S(t,q)
double sq0_ave[3][2]={0},sq1_ave[3][2]={0};//保存了S(0,-q)和S(t,q)的平均值
double phase_real[V]={0},phase_image[V]={0},dssf[2]={0},dssf1[2]={0}; //数组dssf用来保存每一步的S(q,t)*S(-q,0);dssf1用来保存S(q,t)和S(-q,0)的平均值点乘
int cn[4][V]={0};

double ssf_sum[step+2][2]={0};  //在读取每个样本的时候，通过循环对各个样本进行求和，行指标代表的是时间序列/第一列保存了ssf,第二列是ssf*ssf,这里加1是因为有初始值 
double ssf_ave[step+2][2]={0};//ssf_error[step+1]={0};//为了输出方便 
//double spin_sum[step+1][2]={0},spin_ave[step+1]={0},spin_error[step+1]={0};//分别保存了s(t)*s(0)对所有样本求和，对样本求平均，以及误差。
double sqx_sum[step+2][2]={0},sqy_sum[step+2][2]={0},sqz_sum[step+2][2]={0};//分别保存了s(q,t)三个分量的实部和虚部对样本的求和 
//double sqx_ave[step+1][2]={0},sqy_ave[step+1][2]={0},sqz_ave[step+1][2]={0};//分别保存了s(q,t)三个分量的实部和虚部对样本求平均 
double omega=0,ssf_omega[2]={0};
//double drift_sum[step][4]={0};
//矩阵b保存每个时刻各个spin周围的场，根据场和运动方程得到下一刻的自旋构型 
//先固定一个要计算的K点 
table(nn); table1(cn);

double ini_energy=0,energy_t=0;
//double ini_m[3]={0},m_t[3]={0};

//FILE *fp=NULL,*fp1=NULL,*fp2=NULL,*fp3=NULL;
int i=0,j=0,k=0,l=0,s=0;

FILE *fpr[SAMPLE_NUMBER+1]={0}; //定义文件指针，读取各个样本的初始自旋构型 ,注意文件指针是从0开始编号的，自旋样本是从1开始编号的，所以这里有SAMPLE_NUMBER+1 
FILE *fpw[SAMPLE_NUMBER+1]={0}; //定义文件指针，保存各个样本的ssf以及自旋自关联随时间的变化。 
//FILE *fp_drift[SAMPLE_NUMBER+1]={0};   //保存各个样本的守恒量漂移 

FILE *fpa=NULL; //保存dssf对所有样本的平均 
FILE *fp_omega=NULL;
//FILE *fp_da=NULL; //保存RK4更新自旋状态时，平均能量随时间的变化，监测能量漂移，判断计算误差，控制计算steplength

char fileName[MAX_LEN];
char index[MAX_LEN];

char fileName1[MAX_LEN];
char index1[MAX_LEN];

char fileName2[MAX_LEN];
char index2[MAX_LEN];

double wvector[3]={1.0,1.0,1.0};//波矢q=2pi*k(1,1,1)

phase(wvector,phase_real,phase_image,cn); //选择一个固定的K点，首先计算出phase 

for(l=0;l<20;l++)

    {
       printf("%lf         %lf\n",phase_real[l],phase_image[l]);    
    }

for(i=1;i<=SAMPLE_NUMBER;i++)//从样本s1开始依次读取 
{       
    //计算不同的样本前要清零 
    t=0.0;                      
    dssf[0]=0;dssf[1]=0;
    set_zero(sq_0);
    set_zero(sq_1);
    /*
    for(s=0;s<3;s++)
    {
        s0[s]=0;
        ini_m[s]=0;
    }*/     

    //读自旋构型 
    strcpy(fileName,"F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.01_h=0.0\\Ti=0.0\\");
    sprintf(index,"s%d",i);                                    
    strcat(fileName,index);
    strcat(fileName,".txt");
    if((fpr[i]=fopen(fileName,"r"))==NULL)
     {     
          printf("***************************");
          printf("\nCannot open file strike any key exit 1!");

          getch();
          exit(1);
      }

    // if((fp=fopen("F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.1_h=0.0\\Ti=0.0\\s1.txt","r"))==NULL)
    //创建文件夹用于保存每个样本的ssf 
    strcpy(fileName1,"F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.01_h=0.0\\k(1.0_1.0_1.0)\\");
    sprintf(index1,"s%d_ssf",i);
    strcat(fileName1,index1);
    strcat(fileName1,".txt");
    if((fpw[i]=fopen(fileName1,"wt"))==NULL)
     {
          printf("\nCannot open file strike any key exit 2!");
          getch();
          exit(1);
      }
    /*        
    //保存每个样本的能量漂移

    strcpy(fileName2,"F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.02_h=1.0\\k(1.0_1.0_1.0)\\");
    sprintf(index2,"s%d_drift",i);
    strcat(fileName2,index2);
    strcat(fileName2,".txt");
    if((fp_drift[i]=fopen(fileName2,"wt"))==NULL)
     {
          printf("\nCannot open file strike any key exit 3!     %d\n",i);
          getch();
          exit(1);
      }  */    

    for(j=0;j<V;j++)
   {
    fscanf(fpr[i],"%lf",&state[0][j]);
    fscanf(fpr[i],"%lf",&state[1][j]);
    fscanf(fpr[i],"%lf",&state[2][j]);                                
   }

    fclose(fpr[i]); //关闭保存初始自旋构型的文件 

    //计算利用运动方程更新之前体系的初始能量    

     /*ini_energy=tote (state,nn,B); //计算体系的总能量（平均） 
     mom(state,ini_m); //计算体系的总磁矩   */

     //fprintf(fp_drift[i],"%lf          %lf          %lf         %lf         %lf         %lf\n",t,0.0,0.0,0.0,0.0,sqrt(state[0][1024]*state[0][1024]+state[1][1024]*state[1][1024]+state[2][1024]*state[2][1024]));//保存了初始时刻的能量 
    /*
     drift_sum[0][0]=0;
     drift_sum[0][1]=0;
     drift_sum[0][2]=0;
     drift_sum[0][3]=0; */  
    //printf("%lf      %lf      %lf>>>>>\n",state[0][0],state[1][0],state[2][0]);
    //printf("%lf      %lf      %lf>>>>>\n",state[0][2047],state[1][2047],state[2][2047]);

     //首先计算初始时刻的sq 
     cal_sq(state,phase_real,phase_image,sq_1);
     //计算初始时刻的dssf 
     //对于初始时刻,sq_0和sq_1互为复共轭
     Complex_conjugation(sq_0,sq_1);

     cal_dssf(sq_0,sq_1,dssf);
      /*
      printf("**************************\n");
      printf("sq_0:%12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      \n",sq_0[0][0],sq_0[0][1],sq_0[1][0],sq_0[1][1],sq_0[2][0],sq_0[2][1]);
      printf("sq_1:%12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      \n",sq_1[0][0],sq_1[0][1],sq_1[1][0],sq_1[1][1],sq_1[2][0],sq_1[2][1]);
      printf("dssf:%12.10lf      %12.10lf\n",dssf[0],dssf[1]);
      printf("state[][0]:%12.10lf      %12.10lf     %12.10lf\n",state[0][0],state[1][0],state[2][0]);

      printf("ssf_sum[0][0]:%12.10lf      ssf_sum[0][1]:%12.10lf\n",ssf_sum[0][0],ssf_sum[0][1]); 
      printf("**************************\n");        

     for(s=0;s<3;s++)
     {
        s0[s]=state[s][0];
     }
      */
     ssf_sum[0][0]=ssf_sum[0][0]+dssf[0];
     ssf_sum[0][1]=ssf_sum[0][1]+dssf[1];

     //spin_sum[0][0]=spin_sum[0][0]+dot_product(s0,s0);
     //spin_sum[0][1]=spin_sum[0][1]+dot_product(s0,s0)*dot_product(s0,s0);

      sqx_sum[0][0]=sqx_sum[0][0]+sq_1[0][0];
      sqx_sum[0][1]=sqx_sum[0][1]+sq_1[0][1];
      sqy_sum[0][0]=sqy_sum[0][0]+sq_1[1][0];
      sqy_sum[0][1]=sqy_sum[0][1]+sq_1[1][1];
      sqz_sum[0][0]=sqz_sum[0][0]+sq_1[2][0];
      sqz_sum[0][1]=sqz_sum[0][1]+sq_1[2][1];

      //T*t,dot_product(s0,s0)
    fprintf(fpw[i],"%6.4lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf\n",t,dssf[0],dssf[1],sq_1[0][0],sq_1[0][1],sq_1[1][0],sq_1[1][1],sq_1[2][0],sq_1[2][1]);//保存了初始时刻的结构因子 

    for( k=1 ; k < step+1 ; k++)
    {   
      dssf[0]=0;dssf[1]=0;

     /* for(s=0;s<3;s++)
      {
        s1[s]=0;
      } */

      spin_move(state,B);//根据运动方程，利用RK4计算所有spin的状态随时间的演化

      t=t+steplength;
      cal_sq(state,phase_real,phase_image,sq_1);

      //计算新时刻的dssf 
      cal_dssf(sq_0,sq_1,dssf);
      /*
       for(s=0;s<3;s++)
      {
        s1[s]=state[s][0];
      } */

      ssf_sum[k][0]=ssf_sum[k][0]+dssf[0];
      ssf_sum[k][1]=ssf_sum[k][1]+dssf[1];

      //spin_sum[k][0]=spin_sum[k][0]+dot_product(s0,s1);
      //spin_sum[k][1]=spin_sum[k][1]+dot_product(s0,s1)*dot_product(s0,s1);

      sqx_sum[k][0]=sqx_sum[k][0]+sq_1[0][0];
      sqx_sum[k][1]=sqx_sum[k][1]+sq_1[0][1];
      sqy_sum[k][0]=sqy_sum[k][0]+sq_1[1][0];
      sqy_sum[k][1]=sqy_sum[k][1]+sq_1[1][1];
      sqz_sum[k][0]=sqz_sum[k][0]+sq_1[2][0];
      sqz_sum[k][1]=sqz_sum[k][1]+sq_1[2][1];


      /*
      printf("**************************\n");
      printf("sq_0:%12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      \n",sq_0[0][0],sq_0[0][1],sq_0[1][0],sq_0[1][1],sq_0[2][0],sq_0[2][1]);
      printf("sq_1:%12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      %12.10lf      \n",sq_1[0][0],sq_1[0][1],sq_1[1][0],sq_1[1][1],sq_1[2][0],sq_1[2][1]);
      printf("dssf:%12.10lf      %12.10lf\n",dssf[0],dssf[1]);
      printf("state[][0]:%12.10lf      %12.10lf     %12.10lf\n",state[0][0],state[1][0],state[2][0]);

      printf("ssf_sum[0][0]:%12.10lf      ssf_sum[0][1]:%12.10lf\n",ssf_sum[0][0],ssf_sum[0][1]);

      printf("**************************\n");       */

    //fprintf(fpw[i],"%6.4lf          %16.14lf          %16.14lf          %16.14lf          %16.14lf\n",t,dssf[0],dssf[1],T*t,dot_product(s0,s1));//保存了新时刻的结构因子
    fprintf(fpw[i],"%6.4lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf   %16.14lf\n",t,dssf[0],dssf[1],sq_1[0][0],sq_1[0][1],sq_1[1][0],sq_1[1][1],sq_1[2][0],sq_1[2][1]);//保存了初始时刻的结构因子 

    /*     
    m_t[0]=0;m_t[1]=0;m_t[2]=0;

    energy_t=tote (state,nn,B); //计算体系的总能量 

    mom(state,m_t);

    fprintf(fp_drift[i],"%lf           %12.10lf          %12.10lf         %12.10lf          %12.10lf          %12.10lf\n",t,(energy_t-ini_energy),(m_t[0]-ini_m[0]),(m_t[1]-ini_m[1]),(m_t[2]-ini_m[2]),sqrt(state[0][0]*state[0][0]+state[1][0]*state[1][0]+state[2][0]*state[2][0]));//保存了初始时刻的能量 

     drift_sum[k][0]=drift_sum[k][0]+(energy_t-ini_energy);
     drift_sum[k][1]=drift_sum[k][1]+(m_t[0]-ini_m[0]);
     drift_sum[k][2]=drift_sum[k][2]+(m_t[1]-ini_m[1]);
     drift_sum[k][3]=drift_sum[k][3]+(m_t[2]-ini_m[2]); 
    */     
   }

    printf("%d     %d\n",i,fclose(fpw[i])); 
    //fclose(fp_drift[i]);      
}

if((fpa=fopen("F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.01_h=0.0\\k(1.0_1.0_1.0)\\dssf.txt","wt"))==NULL)
{
        printf("\nCannot open file strike any key exit 4!");
        getch();
        exit(1);
}

//首先计算S(-q,0)的平均值   
sq1_ave[0][0]=sqx_sum[0][0]/(SAMPLE_NUMBER);
sq1_ave[0][1]=sqx_sum[0][1]/(SAMPLE_NUMBER);
sq1_ave[1][0]=sqy_sum[0][0]/(SAMPLE_NUMBER);
sq1_ave[1][1]=sqy_sum[0][1]/(SAMPLE_NUMBER);
sq1_ave[2][0]=sqz_sum[0][0]/(SAMPLE_NUMBER);
sq1_ave[2][1]=sqz_sum[0][1]/(SAMPLE_NUMBER);

printf("%lf      %lf     %lf      %lf       %lf       %lf\n",sq1_ave[0][0],sq1_ave[0][1],sq1_ave[1][0],sq1_ave[1][1],sq1_ave[2][0],sq1_ave[2][1]);

Complex_conjugation(sq0_ave,sq1_ave);

ssf_ave[0][0]=1.0*ssf_sum[0][0]/(SAMPLE_NUMBER);
ssf_ave[0][1]=1.0*ssf_sum[0][1]/(SAMPLE_NUMBER);

printf("ssf_ave:%lf      %lf\n",ssf_ave[0][0],ssf_ave[0][1]);

cal_dssf(sq0_ave,sq1_ave,dssf1);

ssf_ave[0][0]=ssf_ave[0][0]-dssf1[0];
ssf_ave[0][1]=ssf_ave[0][1]-dssf1[1];


fprintf(fpa,"%6.4lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf\n",0.0,ssf_ave[0][0],ssf_ave[0][1],sq1_ave[0][0],sq1_ave[0][1],sq1_ave[1][0],sq1_ave[1][1],sq1_ave[2][0],sq1_ave[2][1]);


for(i=1;i<step+1;i++)
{
    dssf1[0]=0;dssf1[1]=0;
    set_zero(sq1_ave);
    ssf_ave[i][0]=1.0*ssf_sum[i][0]/(SAMPLE_NUMBER);
    ssf_ave[i][1]=1.0*ssf_sum[i][1]/(SAMPLE_NUMBER);

    if(i==0)
    {

    printf("%lf      %lf     %lf      %lf       %lf       %lf\n",sq1_ave[0][0],sq1_ave[0][1],sq1_ave[1][0],sq1_ave[1][1],sq1_ave[2][0],sq1_ave[2][1]);

    }
    //ssf_error[i]=1.0*sqrt((ssf_sum[i][1]-1.0*(SAMPLE_NUMBER)*ssf_ave[i]*ssf_ave[i])/(SAMPLE_NUMBER*(SAMPLE_NUMBER-1.0)));

    //spin_ave[i]=1.0*spin_sum[i][0]/(SAMPLE_NUMBER);
    //spin_error[i]=1.0*sqrt((spin_sum[i][1]-1.0*(SAMPLE_NUMBER)*spin_ave[i]*spin_ave[i])/(SAMPLE_NUMBER*(SAMPLE_NUMBER-1.0)));

    sq1_ave[0][0]=sqx_sum[i][0]/(SAMPLE_NUMBER);
    sq1_ave[0][1]=sqx_sum[i][1]/(SAMPLE_NUMBER);
    sq1_ave[1][0]=sqy_sum[i][0]/(SAMPLE_NUMBER);
    sq1_ave[1][1]=sqy_sum[i][1]/(SAMPLE_NUMBER);
    sq1_ave[2][0]=sqz_sum[i][0]/(SAMPLE_NUMBER);
    sq1_ave[2][1]=sqz_sum[i][1]/(SAMPLE_NUMBER);

    //这里要检查S(q,t)确实是时间平移不变的?????好像并不是 

    //计算<S(q,t)S(-q,0)>-<S(q,t)><S(-q,0)> 
    cal_dssf(sq0_ave,sq1_ave,dssf1);


    ssf_ave[i][0]=ssf_ave[i][0]-dssf1[0];
    ssf_ave[i][1]=ssf_ave[i][1]-dssf1[1];

    fprintf(fpa,"%6.4lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf   %16.12lf\n",steplength*i,ssf_ave[i][0],ssf_ave[i][1],sq1_ave[0][0],sq1_ave[0][1],sq1_ave[1][0],sq1_ave[1][1],sq1_ave[2][0],sq1_ave[2][1]);

}

printf("fclose(fpa):%d\n",fclose(fpa));

if((fp_omega=fopen("F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.01_h=0.0\\k(1.0_1.0_1.0)\\dssf_omega.txt","wt"))==NULL)
{
        printf("\nCannot open file strike any key exit 5!");
        getch();
        exit(1);
}

for(omega=-4.0;omega<=4.01;)
{
    ssf_omega[0]=0;ssf_omega[1]=0;

    for(i=0;i<step+1;i++)
    {
        ssf_omega[0]=ssf_omega[0]+ssf_ave[i][0]*cos(omega*i*steplength)*steplength-ssf_ave[i][1]*sin(omega*i*steplength)*steplength;
        ssf_omega[1]=ssf_omega[1]+ssf_ave[i][0]*sin(omega*i*steplength)*steplength+ssf_ave[i][1]*cos(omega*i*steplength)*steplength;

    }

    fprintf(fp_omega,"%lf         %lf         %lf          %lf       %lf\n",omega,2*ssf_omega[0],ssf_omega[1],(omega/(4*PI*PI*dot_product(wvector,wvector))),2*(4*PI*PI*ssf_omega[0]*dot_product(wvector,wvector)/T));

    omega=omega+0.005;
}

printf("fclose(fp_omega):%d\n",fclose(fp_omega));   
printf("%lf\n",dot_product(wvector,wvector));

/*最后一段是用来计算守恒量的漂移随时间的变化*/  

/*      
   if((fp_da=fopen("F:\\Pyrochlore_Data\\pyrochlore_smaple\\Tem=0.02_h=1.0\\k(1.0_1.0_1.0)\\drift_sl=0.025.txt","wt"))==NULL)
    {
          printf("\nCannot open file strike any key exit 6!");
          getch();
          exit(1);
    }

   for( k=0 ; k < step+0.5 ; k++)
  {
    fprintf(fp_da,"%lf         %12.10lf           %12.10lf         %12.10lf           %12.10lf\n",steplength*k,(drift_sum[k][0]/(SAMPLE_NUMBER)),(drift_sum[k][1]/(SAMPLE_NUMBER)),(drift_sum[k][2]/(SAMPLE_NUMBER)),(drift_sum[k][3]/(SAMPLE_NUMBER)));
 }

  fclose(fp_da);
  */                  
  return 0;

}

也怀疑是自己写的循环有问题。但是单独拿出来运行又正常了。不明白是怎么回事。

  • 写回答

0条回答 默认 最新

    报告相同问题？

    悬赏问题

    • ¥15 乘性高斯噪声在深度学习网络中的应用
    • ¥15 运筹学排序问题中的在线排序
    • ¥15 关于docker部署flink集成hadoop的yarn，请教个问题 flink启动yarn-session.sh连不上hadoop，这个整了好几天一直不行，求帮忙看一下怎么解决
    • ¥30 求一段fortran代码用IVF编译运行的结果
    • ¥15 深度学习根据CNN网络模型，搭建BP模型并训练MNIST数据集
    • ¥15 C++ 头文件/宏冲突问题解决
    • ¥15 用comsol模拟大气湍流通过底部加热（温度不同）的腔体
    • ¥50 安卓adb backup备份子用户应用数据失败
    • ¥20 有人能用聚类分析帮我分析一下文本内容嘛
    • ¥30 python代码，帮调试，帮帮忙吧