gnmlast.c

gauss network model方法研究蛋白质折叠的程序

#include <stdio.h>
#include <math.h>
#include <string.h>

main()
{
  char filein[40],fileout[40],filetwo[40];
  FILE *fpi,*fpo,*fpt;
  float coor[300][3];  //put the CA coor of the protein
  float kirchhoff[300][300];  //is the kirchhoff matix
  int i,j,k,l,m,presnum;
  int resnum[300];
  char line[80];
  char test[5],test1[3];
  float max,dist;
  int maxi,maxj;
   float dctemp[300][300];
   float cmatix[300][300];
   float benzhengxl[300][300];
   float benzhengxltemp[300][300];
   float epslo;
   float phi;
   for(i=0;i<300;i++)
     {
       for(j=0;j<300;j++)
         {benzhengxl[i][j]=0.0;
          benzhengxltemp[i][j]=0.0;}
       benzhengxltemp[i][i]=1.0;
     }   
//the following is read the CA coor of the protein
  scanf("%s\n",filein);
  printf("%s\n",filein);
  scanf("%s\n",fileout);
  printf("%s\n",fileout);
  scanf("%s\n",filetwo);
  printf("%s\n",filetwo);

  if((fpi=fopen(filein,"r"))==NULL)
    {printf("can not open the input file");
     exit(0);
    }
   if((fpo=fopen(fileout,"w"))==NULL)
    {printf("can not open the output file");
     exit(0);
    }
    if((fpt=fopen(filetwo,"w"))==NULL)
    {printf("can not open the output file");
     exit(0);
    }
  fgets(line,80,fpi);
presnum=0;
  while(!feof(fpi))
   {//printf("ok!\n");
     strncpy(test,line,4);test[4]='\0';
     strncpy(test1,line+13,2);test1[2]='\0';
     if((strcmp(test,"ATOM")==0&&(strcmp(test1,"CA")==0)))
       { // printf("ok\n");
          sscanf(line+28,"%f",&coor[presnum][0]);
          sscanf(line+38,"%f",&coor[presnum][1]);
          sscanf(line+46,"%f",&coor[presnum][2]);
          sscanf(line+22,"%d",&resnum[presnum]);
          presnum++;
        }
      fgets(line,80,fpi);
    }

//the following is calculating the Kirchhoff matix
  presnum=presnum-1;
printf("presnum=%d\n",presnum);
  for(i=0;i<presnum;i++)
    {
      for(j=i+1;j<=presnum;j++)
         {
           dist=sqrt((coor[j][0]-coor[i][0])*(coor[j][0]-coor[i][0])+(coor[j][
1]-coor[i][1])*(coor[j][1]-coor[i][1])+(coor[j][2]-coor[i][2])*(coor[j][2]-coor[i][2]));
          // printf("dist=%f\n",dist);
            if(dist<=7.3)
              {kirchhoff[i][j]=-1.0;}
            else
              {kirchhoff[i][j]=0.0;}
            kirchhoff[j][i]=kirchhoff[i][j];
           }
     }
 
  for(i=0;i<=presnum;i++)
    {
      kirchhoff[i][i]=0;
      for(j=0;j<=presnum;j++)
       {
         if(j!=i)
           {kirchhoff[i][i]=kirchhoff[i][i]-kirchhoff[i][j];}
        }
     }
   for(i=0;i<=presnum;i++)
   { 
    printf("kirchhoff[%d][%d]=%f\n",i,i,kirchhoff[i][i]);
   }
/*kirchhoff[0][0]=2.0;kirchhoff[0][1]=-1.0;kirchhoff[0][2]=0.0;kirchhoff[1][0]=-1.0;kirchhoff[1][1]=2.0;kirchhoff[1][2]=-1.0;
kirchhoff[2][0]=0.0;kirchhoff[2][1]=-1.0;kirchhoff[2][2]=2.0;presnum=2;*/
//the following is cal the eigntors and eignvectors of the matix

 loop:epslo=0.0;
   maxi=0;maxj=1;
   max=kirchhoff[0][1];
   for(i=0;i<=presnum;i++)
    {
       for(j=0;j<=presnum;j++)
         {
           if(j!=i)
             {
               if(fabs(kirchhoff[i][j])>fabs(max))
                 {maxi=i;maxj=j;max=kirchhoff[i][j];}
             }
         }
     }
    printf("i=%d,j=%d\n",maxi,maxj);
    phi=atan(2*kirchhoff[maxi][maxj]/(kirchhoff[maxi][maxi]-kirchhoff[maxj][maxj]
))/2;
    dctemp[maxi][maxi]=kirchhoff[maxi][maxi]*cos(phi)*cos(phi)+kirchhoff[maxj][maxj]*sin(phi)*sin(phi)+2*kirchhoff[maxi][maxj]*cos(phi)*sin(phi);

    dctemp[maxj][maxj]=kirchhoff[maxi][maxi]*sin(phi)*sin(phi)+kirchhoff[maxj][maxj]*cos(phi)*cos(phi)-2*kirchhoff[maxi][maxj]*cos(phi)*sin(phi);
    for(l=0;l<=presnum;l++)
      {if((l!=maxi)&&(l!=maxj))
       dctemp[maxi][l]=kirchhoff[maxi][l]*cos(phi)+kirchhoff[maxj][l]*sin(phi);
       dctemp[l][maxi]=dctemp[maxi][l];
      }
     for(l=0;l<=presnum;l++)
      {if((l!=maxi)&&(l!=maxj))
       dctemp[maxj][l]=kirchhoff[maxj][l]*cos(phi)-kirchhoff[maxi][l]*sin(phi);
       dctemp[l][maxj]=dctemp[maxj][l];
      }
       for(l=0;l<=presnum;l++)
        {
           for(m=0;m<=presnum;m++)
             {if((l!=maxi)&&(l!=maxj)&&(m!=maxi)&&(m!=maxj))
               {dctemp[l][m]=kirchhoff[l][m];
                dctemp[m][l]=dctemp[l][m];}
              }
         }
      dctemp[maxi][maxj]=(kirchhoff[maxj][maxj]-kirchhoff[maxi][maxi])*sin(2*phi)/2+kirchhoff[maxi][maxj]*(cos(phi)*cos(phi)-sin(phi)*sin(phi));
      dctemp[maxj][maxi]=dctemp[maxi][maxj];
      for(i=0;i<=presnum;i++)
       {
        for(j=0;j<=presnum;j++)
         {kirchhoff[i][j]=dctemp[i][j];}
        }

 //the following is cal the benzhengxiangliang matix
      for(i=0;i<=presnum;i++)
         {
           for(j=0;j<=presnum;j++)
              {
                if(j!=i)
                 {cmatix[i][j]=0.0;}
              }
            cmatix[i][i]=1.0;
          }
       cmatix[maxi][maxi]=cos(phi);
       cmatix[maxj][maxj]=cos(phi);
       cmatix[maxi][maxj]=-sin(phi);
       cmatix[maxj][maxi]=sin(phi);
      for(i=0;i<=presnum;i++)
        {
           for(j=0;j<=presnum;j++)
             {benzhengxl[i][j]=0.0;}
         }
       for(i=0;i<=presnum;i++)
         {
            for(j=0;j<=presnum;j++)
              {
                 for(k=0;k<=presnum;k++)
                   {benzhengxl[i][j]=benzhengxl[i][j]+benzhengxltemp[i][k]*cmatix[k][j];}
            // printf("benzhengxl[%d][%d]=%f\n",i,j,benzhengxl[i][j]);
              }
          }
          for(i=0;i<=presnum;i++)
        {
            for(j=0;j<=presnum;j++)
                 {benzhengxltemp[i][j]=benzhengxl[i][j];}
         }
       for(i=0;i<=presnum;i++)
          {
             for(j=0;j<=presnum;j++)
               {
                  if (j!=i)
                   {epslo=epslo+dctemp[i][j]*dctemp[i][j];}
                }
           }
         if(epslo>0.0000000001)
          {goto loop;}
         for(i=0;i<=presnum;i++)
            { fprintf(fpo,"%5d%12.6f\n",i,dctemp[i][i]);
              printf("dctemp[%d][%d]=%f\n",i,i,dctemp[i][i]);}
         for(i=0;i<=presnum;i++)
            {
               for(j=0;j<=presnum;j++)
                  {fprintf(fpt,"%5d%5d%12.6f\n",i,j,benzhengxl[i][j]);
                   printf("benzhengxl[%d][%d]=%f\n",i,j,benzhengxl[i][j]); }

            }
          fclose(fpi);
          fclose(fpo);
          fclose(fpt);
}