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

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


# define  N   8
# define  pi  3.14159265

double th0=1.57;
double ph0=2.0944;//0.873;
double l=4.0;

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

int cart2dof();

int main()
{
  int i, j, k;  
  double dx, dy, dz, dr;
  double A, B, C, D, E, F, G, AA, BB, CC, DD, EE, FF, GG, CV1, CV2, CV3, CV4;
  double xm, ym, zm;
  FILE  *fp;

  ym=-cos(ph0)*sin(th0)*cos(th0/2)+sin(th0/2)-cos(th0)*sin(th0/2);
  ym=ym*l;

  xm=(l*cos(ph0)*sin(th0)+ym*cos(th0/2))/sin(th0/2);
  xm=-xm;

  zm=l*l-xm*xm-l*l*cos(th0/2)*cos(th0/2)-2*xm*l*cos(th0/2)-ym*ym-l*l*sin(th0/2)*sin(th0/2)+2*ym*l*sin(th0/2);
  zm=sqrt(zm);

  A=xm/2+l*cos(th0/2);
  B=ym/2-l*sin(th0/2);
  C=zm/2;
  
  E=-2*B*xm/(sqrt(A*A+B*B+C*C))+3*l*sin(th0/2);
  F=-2*C*xm/(sqrt(A*A+B*B+C*C));
  G=-2*A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2);

  AA=-A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2)-xm;
  BB=-B*xm/(sqrt(A*A+B*B+C*C))+2*l*sin(th0/2)-ym;
  CC=-C*xm/(sqrt(A*A+B*B+C*C))-zm;



  dr=2*sqrt(ym*ym+zm*zm);

  printf("dr: %g\n",dr);

  fp=fopen("distance","a");
  fprintf(fp,"%g  %g  %g\n",th0,ph0,dr);
  fclose(fp);

  x[6]=-2*AA*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2xm;y[6]=-2*BB*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2ym;z[6]=-2*CC*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2zm;
  x[5]=G;y[5]=E;z[5]=F;
  x[4]=xm; y[4]=ym; z[4]=zm;
  x[3]=-l*cos(th0/2); y[3]=l*sin(th0/2); z[3]=0;
  x[2]=0; y[2]=0; z[2]=0;
  x[1]=-l*cos(th0/2); y[1]=-l*sin(th0/2); z[1]=0;
  x[0]=xm; y[0]=-ym; z[0]=-zm;

  DD=((x[6]+x[4])/2-x[5])*((x[6]+x[4])/2-x[5]);
  EE=((y[6]+y[4])/2-y[5])*((y[6]+y[4])/2-y[5]);
  FF=((z[6]+z[4])/2-z[5])*((z[6]+z[4])/2-z[5]);

  GG=sqrt(DD+EE+FF);

  x[7]=(x[6]+x[4]-2*x[5])*(abs(xm))/GG+2*x[5]-x[3];

  for(i=0;i<N-1;i++)
  {
     j=i+1;
     dx=x[j]-x[i];
     dy=y[j]-y[i];
     dz=z[j]-z[i];
     dr=sqrt(dx*dx+dy*dy+dz*dz);
     printf("bond length  i=%d  %g\n",i+1,dr);
  }

  cart2dof();

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

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

  fp=fopen("conf","w");
  for(i=0;i<N;i++)
  {
     fprintf(fp,"%.5f  %.5f  %.5f\n",x[i],y[i],z[i]);
  }
  fclose(fp);


  CV3=x[5]-xm;//(H5,CV3>0)
  
  
  CV4=x[7]-xm;//(H7,CV4<0)

  fp=fopen("r57","a");
  fprintf(fp,"CV1 %g\n",CV3);
  fprintf(fp,"CV2 %g\n",CV4);

  /*CV1=G;  // H5 (G<0)
  CV2=0.5*cos(th0)/(1-cos(th0/2)*cos(th0/2))-cos(ph0);  // H3, CV2 <0 

  fp=fopen("r35","a");
  fprintf(fp,"CV1 %g\n",CV1);
  fprintf(fp,"CV2 %g\n",CV2);
  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<=-0.999999) {tmp1=-0.999999; ok=0;}
    if (tmp1>=0.999999) {tmp1=0.999999; ok=0;}
    if (acos(tmp1)==0.0) {ok=0;}
    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<=-0.999999) {tmp1=-0.999999; ok=0;}
     if(tmp1>=0.999999)  {tmp1=0.999999; ok=0;}
      
     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;
}

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


# define  N   8
# define  pi  3.14159265

double th0=1.57;
double ph0=2.0944;//0.873;
double l=4.0;

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

int cart2dof();

int main()
{
  int i, j, k;  
  double dx, dy, dz, dr;
  double A, B, C, D, E, F, G, AA, BB, CC, DD, EE, FF, GG, CV1, CV2, CV3, CV4;
  double xm, ym, zm;
  FILE  *fp;

  ym=-cos(ph0)*sin(th0)*cos(th0/2)+sin(th0/2)-cos(th0)*sin(th0/2);
  ym=ym*l;

  xm=(l*cos(ph0)*sin(th0)+ym*cos(th0/2))/sin(th0/2);
  xm=-xm;

  zm=l*l-xm*xm-l*l*cos(th0/2)*cos(th0/2)-2*xm*l*cos(th0/2)-ym*ym-l*l*sin(th0/2)*sin(th0/2)+2*ym*l*sin(th0/2);
  zm=sqrt(zm);

  A=xm/2+l*cos(th0/2);
  B=ym/2-l*sin(th0/2);
  C=zm/2;
  
  E=-2*B*xm/(sqrt(A*A+B*B+C*C))+3*l*sin(th0/2);
  F=-2*C*xm/(sqrt(A*A+B*B+C*C));
  G=-2*A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2);

  AA=-A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2)-xm;
  BB=-B*xm/(sqrt(A*A+B*B+C*C))+2*l*sin(th0/2)-ym;
  CC=-C*xm/(sqrt(A*A+B*B+C*C))-zm;



  dr=2*sqrt(ym*ym+zm*zm);

  printf("dr: %g\n",dr);

  fp=fopen("distance","a");
  fprintf(fp,"%g  %g  %g\n",th0,ph0,dr);
  fclose(fp);

  x[6]=-2*AA*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2xm;y[6]=-2*BB*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2ym;z[6]=-2*CC*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2zm;
  x[5]=G;y[5]=E;z[5]=F;
  x[4]=xm; y[4]=ym; z[4]=zm;
  x[3]=-l*cos(th0/2); y[3]=l*sin(th0/2); z[3]=0;
  x[2]=0; y[2]=0; z[2]=0;
  x[1]=-l*cos(th0/2); y[1]=-l*sin(th0/2); z[1]=0;
  x[0]=xm; y[0]=-ym; z[0]=-zm;

  DD=((x[6]+x[4])/2-x[5])*((x[6]+x[4])/2-x[5]);
  EE=((y[6]+y[4])/2-y[5])*((y[6]+y[4])/2-y[5]);
  FF=((z[6]+z[4])/2-z[5])*((z[6]+z[4])/2-z[5]);

  GG=sqrt(DD+EE+FF);

  x[7]=(x[6]+x[4]-2*x[5])*(abs(xm))/GG+2*x[5]-x[3];

  for(i=0;i<N-1;i++)
  {
     j=i+1;
     dx=x[j]-x[i];
     dy=y[j]-y[i];
     dz=z[j]-z[i];
     dr=sqrt(dx*dx+dy*dy+dz*dz);
     printf("bond length  i=%d  %g\n",i+1,dr);
  }

  cart2dof();

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

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

  fp=fopen("conf","w");
  for(i=0;i<N;i++)
  {
     fprintf(fp,"%.5f  %.5f  %.5f\n",x[i],y[i],z[i]);
  }
  fclose(fp);


  CV3=x[5]-xm;//(H5,CV3>0)
  
  
  CV4=x[7]-xm;//(H7,CV4<0)

  fp=fopen("r57","a");
  fprintf(fp,"CV1 %g\n",CV3);
  fprintf(fp,"CV2 %g\n",CV4);

  /*CV1=G;  // H5 (G<0)
  CV2=0.5*cos(th0)/(1-cos(th0/2)*cos(th0/2))-cos(ph0);  // H3, CV2 <0 

  fp=fopen("r35","a");
  fprintf(fp,"CV1 %g\n",CV1);
  fprintf(fp,"CV2 %g\n",CV2);
  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<=-0.999999) {tmp1=-0.999999; ok=0;}
    if (tmp1>=0.999999) {tmp1=0.999999; ok=0;}
    if (acos(tmp1)==0.0) {ok=0;}
    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<=-0.999999) {tmp1=-0.999999; ok=0;}
     if(tmp1>=0.999999)  {tmp1=0.999999; ok=0;}
      
     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;
}

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


# define  N   8
# define  pi  3.14159265

double th0=1.57;
double ph0=2.0944;//0.873;
double l=4.0;

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

int cart2dof();

int main()
{
  int i, j, k;  
  double dx, dy, dz, dr;
  double A, B, C, D, E, F, G, AA, BB, CC, DD, EE, FF, GG, CV1, CV2, CV3, CV4;
  double xm, ym, zm;
  FILE  *fp;

  ym=-cos(ph0)*sin(th0)*cos(th0/2)+sin(th0/2)-cos(th0)*sin(th0/2);
  ym=ym*l;

  xm=(l*cos(ph0)*sin(th0)+ym*cos(th0/2))/sin(th0/2);
  xm=-xm;

  zm=l*l-xm*xm-l*l*cos(th0/2)*cos(th0/2)-2*xm*l*cos(th0/2)-ym*ym-l*l*sin(th0/2)*sin(th0/2)+2*ym*l*sin(th0/2);
  zm=sqrt(zm);

  A=xm/2+l*cos(th0/2);
  B=ym/2-l*sin(th0/2);
  C=zm/2;
  
  E=-2*B*xm/(sqrt(A*A+B*B+C*C))+3*l*sin(th0/2);
  F=-2*C*xm/(sqrt(A*A+B*B+C*C));
  G=-2*A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2);

  AA=-A*xm/(sqrt(A*A+B*B+C*C))-l*cos(th0/2)-xm;
  BB=-B*xm/(sqrt(A*A+B*B+C*C))+2*l*sin(th0/2)-ym;
  CC=-C*xm/(sqrt(A*A+B*B+C*C))-zm;



  dr=2*sqrt(ym*ym+zm*zm);

  printf("dr: %g\n",dr);

  fp=fopen("distance","a");
  fprintf(fp,"%g  %g  %g\n",th0,ph0,dr);
  fclose(fp);

  x[6]=-2*AA*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2xm;y[6]=-2*BB*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2ym;z[6]=-2*CC*xm/(sqrt(AA*AA+BB*BB+CC*CC))+2zm;
  x[5]=G;y[5]=E;z[5]=F;
  x[4]=xm; y[4]=ym; z[4]=zm;
  x[3]=-l*cos(th0/2); y[3]=l*sin(th0/2); z[3]=0;
  x[2]=0; y[2]=0; z[2]=0;
  x[1]=-l*cos(th0/2); y[1]=-l*sin(th0/2); z[1]=0;
  x[0]=xm; y[0]=-ym; z[0]=-zm;

  DD=((x[6]+x[4])/2-x[5])*((x[6]+x[4])/2-x[5]);
  EE=((y[6]+y[4])/2-y[5])*((y[6]+y[4])/2-y[5]);
  FF=((z[6]+z[4])/2-z[5])*((z[6]+z[4])/2-z[5]);

  GG=sqrt(DD+EE+FF);

  x[7]=(x[6]+x[4]-2*x[5])*(abs(xm))/GG+2*x[5]-x[3];

  for(i=0;i<N-1;i++)
  {
     j=i+1;
     dx=x[j]-x[i];
     dy=y[j]-y[i];
     dz=z[j]-z[i];
     dr=sqrt(dx*dx+dy*dy+dz*dz);
     printf("bond length  i=%d  %g\n",i+1,dr);
  }

  cart2dof();

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

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

  fp=fopen("conf","w");
  for(i=0;i<N;i++)
  {
     fprintf(fp,"%.5f  %.5f  %.5f\n",x[i],y[i],z[i]);
  }
  fclose(fp);


  CV3=x[5]-xm;//(H5,CV3>0)
  
  
  CV4=x[7]-xm;//(H7,CV4<0)

  fp=fopen("r57","a");
  fprintf(fp,"CV1 %g\n",CV3);
  fprintf(fp,"CV2 %g\n",CV4);

  /*CV1=G;  // H5 (G<0)
  CV2=0.5*cos(th0)/(1-cos(th0/2)*cos(th0/2))-cos(ph0);  // H3, CV2 <0 

  fp=fopen("r35","a");
  fprintf(fp,"CV1 %g\n",CV1);
  fprintf(fp,"CV2 %g\n",CV2);
  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<=-0.999999) {tmp1=-0.999999; ok=0;}
    if (tmp1>=0.999999) {tmp1=0.999999; ok=0;}
    if (acos(tmp1)==0.0) {ok=0;}
    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<=-0.999999) {tmp1=-0.999999; ok=0;}
     if(tmp1>=0.999999)  {tmp1=0.999999; ok=0;}
      
     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;
}

  • 6
    点赞
  • 9
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值