C语言——一些处理脚本(6)——(待完善)

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

# define N   30

# define max(A,B) ((A)>(B)?(A):(B)) 
# define min(A,B) ((A)<(B)?(A):(B))

double kth=10;
double kph1=5;

double th0=1.570;//1.832;
double ph0=0.873;


# define L0   14
# define L1   15

double cthmin, sthmin;
int therr=0;

double pi; 
double x[N], y[N], z[N];
double b[N];                       /* pseudo bond lenghts (0...N-1)                   */
double th[N];                      /* pseudo bond angles (1...N-2)                    */
double bx[N],by[N],bz[N];          /* pseudo bond vectors                             */
double ph[N];                      /* pseudo torsion angles (1...N-3)                 */
double sx[N],sy[N],sz[N];          /* auxilary vector                                 */
double Eben[N], Etor[N];

int cart2dof();
double bend();
double torsion();

int main(void)
{
  int i, j, k;  
  double Eth=0, Etau=0;
  char Istring[50];
  FILE *fp, *fp1;

   pi=acos(-1.);
   cthmin=cos(pi/180.);
   sthmin=sin(pi/180.);

  printf("k=");
  scanf("%d",&k);
  sprintf(Istring,"snp%d",k);
    
  fp=fopen(Istring,"r");
  for(i=0;i<N;i++)
  {
      fscanf(fp,"%lf %lf %lf",&x[i],&y[i],&z[i]);
  }
  fclose(fp);
  
  cart2dof();
  Eth=bend();
  Etau=torsion();

  fp=fopen("90theta","w");
  fprintf(fp,"%g\n\n",Eth);
  for(i=1;i<N-1;i++)
  {
     fprintf(fp,"%d  %g  %g  %g\n",i+1,th[i],Eben[i],th0);
  }  
  fclose(fp);

  fp=fopen("90tau","w");
  fprintf(fp,"%g\n\n",Etau);
  for(i=1;i<N-2;i++)
  {
     fprintf(fp,"%d  %g  %g  %g\n",i+1,ph[i],Etor[i],ph0);
  }  
  fclose(fp);


  return 0;
}

int cart2dof()   //Caculation of bond, bond angles, torsion angles.
{
  int i,j,ok=1;
  double b1x,b1y,b1z,b1;
  double b2x,b2y,b2z;
  double ux,uy,uz,u;
  double tmp1;

  /* bond vectors */
  for(i=0;i<N-1;i++) 
  {
    j=i+1;
    b1x=x[j]-x[i];
    b1y=y[j]-y[i];
    b1z=z[j]-z[i];
    b1=sqrt(b1x*b1x+b1y*b1y+b1z*b1z);
    
    bx[i]=b1x/b1;
    by[i]=b1y/b1;
    bz[i]=b1z/b1;
    b[i]=b1;
    if(b1>5.0 || b1<3.0)
    {
       ok=0; 
    }
  }    

  /* bond angles */
  for(i=1;i<N-1;i++) 
  {
    b1x=bx[i-1];
    b1y=by[i-1];
    b1z=bz[i-1];

    b2x=bx[i];
    b2y=by[i];
    b2z=bz[i];
    
    tmp1=b1x*b2x+b1y*b2y+b1z*b2z;
    if (tmp1<=-1.0) {tmp1=-0.999999; ok=0;}// printf("cart2dof 2: tmp1 %f\n",tmp1);}
    if (tmp1>=1.0) {tmp1=0.999999; ok=0;}// printf("cart2dof 3: tmp1 %f\n",tmp1);}
    if (acos(tmp1)==0.0) {ok=0;}// printf("cart2dof 4: tmp1 %f\n",tmp1);}
    th[i]=pi-acos(tmp1);

    ux=b1y*b2z-b1z*b2y;
    uy=b1z*b2x-b1x*b2z;
    uz=b1x*b2y-b1y*b2x;
    u=sqrt(ux*ux+uy*uy+uz*uz);

    /* auxiliary vectors */ 
    sx[i]=ux/u;
    sy[i]=uy/u;
    sz[i]=uz/u;
  }
 
  /* torsion angles */
  for(i=1;i<N-2;i++) 
  {
    j=i+1;
    
    tmp1=sx[i]*sx[j]+sy[i]*sy[j]+sz[i]*sz[j];
    if(tmp1<=-1.0) {tmp1=-0.999999; ok=0;}// printf("cart2dof 5: tmp1 %f\n",tmp1);}
    if(tmp1>=1.0)  {tmp1=0.999999; ok=0;}// printf("cart2dof 6: tmp1 %f\n",tmp1);}
      
    if((sx[i]*bx[j]+sy[i]*by[j]+sz[i]*bz[j])> 0)
    {
       ph[i]=acos(tmp1);
    }
    else ph[i]=-acos(tmp1);  
    
  }
  return ok;
}

/*double bend() 
{
  int i,j,k;
  double b1x,b1y,b1z,b1;
  double b2x,b2y,b2z,b2;
  double dix,diy,diz;
  double dkx,dky,dkz;
  double cth,sth,dth;
  double e=0,fben;

  for(i=0;i<N-2;i++) 
  {
    j=i+1;
    k=i+2;
    dth=th[j]-th0;  // alpha helix
    cth=max(cos(th[j]),-cthmin);
    sth=max(sin(th[j]),sthmin);
    if(sin(th[j])<sthmin) therr++;
    
    b1x=bx[i];
    b1y=by[i];
    b1z=bz[i];
    b1=b[i];

    b2x=bx[j];
    b2y=by[j];
    b2z=bz[j];
    b2=b[j];

    e+=dth*dth;
    Eben[j]=kth*dth*dth;
  }

  return kth*e;
}*/

double bend() 
{
  int i,j,k;
  double b1x,b1y,b1z,b1;
  double b2x,b2y,b2z,b2;
  double dix,diy,diz;
  double dkx,dky,dkz;
  double cth,sth,dth;
  double e=0,fben;

  for(i=0;i<L0-2;i++) 
  {
    j=i+1;
    k=i+2;    
    dth=th[j]-th0;  // alpha helix
    cth=max(cos(th[j]),-cthmin);
    sth=max(sin(th[j]),sthmin);
    if(sin(th[j])<sthmin) therr++;
    
    b1x=bx[i];
    b1y=by[i];
    b1z=bz[i];
    b1=b[i];

    b2x=bx[j];
    b2y=by[j];
    b2z=bz[j];
    b2=b[j];

    e+=dth*dth;
    Eben[j]=kth*dth*dth;
  }

  for(i=L1+1;i<N-2;i++) 
  {
    j=i+1;
    k=i+2;    
    dth=th[j]-th0;  // alpha helix
    cth=max(cos(th[j]),-cthmin);
    sth=max(sin(th[j]),sthmin);
    if(sin(th[j])<sthmin) therr++;
    
    b1x=bx[i];
    b1y=by[i];
    b1z=bz[i];
    b1=b[i];

    b2x=bx[j];
    b2y=by[j];
    b2z=bz[j];
    b2=b[j];

    e+=dth*dth;
    Eben[j]=kth*dth*dth;
   }

  return kth*e;
}

/*
double torsion() 
{
  int i,j,k,l;
  double e1=0,e2=0,fph=0;
  double b1,b2,b3;
  double dix,diy,diz;
  double djx,djy,djz;
  double dkx,dky,dkz;
  double dlx,dly,dlz;
  double cth1,cth2,sth1,sth2,dph;

  for(i=0;i<N-3;i++) 
  {
    j=i+1;
    k=i+2;
    l=i+3;

    dph=ph[j]-ph0;
    cth1=max(cos(th[j]),-cthmin);
    sth1=max(sin(th[j]),sthmin);
    cth2=max(cos(th[k]),-cthmin);
    sth2=max(sin(th[k]),sthmin);
    
    b1=b[i]; 
    b2=b[j]; 
    b3=b[k];

    e1+=1-cos(dph);
    Etor[j]=kph1*(1-cos(dph));
    
  }
  return kph1*e1;//+kph2*e2;
}*/


double torsion() //
{
  int i,j,k,l;
  double e1=0,e2=0,fph=0;
  double b1,b2,b3;
  double dix,diy,diz;
  double djx,djy,djz;
  double dkx,dky,dkz;
  double dlx,dly,dlz;
  double cth1,cth2,sth1,sth2,dph;

  for(i=0;i<L0-3;i++) 
  {
    j=i+1;
    k=i+2;
    l=i+3;

    dph=ph[j]-ph0;
    cth1=max(cos(th[j]),-cthmin);
    sth1=max(sin(th[j]),sthmin);
    cth2=max(cos(th[k]),-cthmin);
    sth2=max(sin(th[k]),sthmin);
    
    b1=b[i]; 
    b2=b[j]; 
    b3=b[k];

    e1+=1-cos(dph);
    Etor[j]=kph1*(1-cos(dph));
  }

  for(i=L1+1;i<N-3;i++) 
  {
    j=i+1;
    k=i+2;
    l=i+3;

    dph=ph[j]-ph0;
    cth1=max(cos(th[j]),-cthmin);
    sth1=max(sin(th[j]),sthmin);
    cth2=max(cos(th[k]),-cthmin);
    sth2=max(sin(th[k]),sthmin);
    
    b1=b[i]; 
    b2=b[j]; 
    b3=b[k];

    e1+=1-cos(dph);
    Etor[j]=kph1*(1-cos(dph)); 
  }
  return kph1*e1;
}