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

#define NATOMX 160
#define NATOMY 300



/* lenght in picometer !!!! if not, round-off errors !!!! */

void main(int argc, char *argv[]){

float atoms[NATOMX*NATOMY][2],rot[NATOMX*NATOMY][2],rt1[NATOMX*NATOMY][3],rt2[NATOMX*NATOMY][3];
float sheet[NATOMX*NATOMY][2],a=142,xmin,xmax,ymin,ymax,cir,sig=340.0,eps=0.012,evdw,energy,potcut,cut=3.2,rad;
float phi,teta, alpha,r,tube1[NATOMX*NATOMY][3],delta,tube2[NATOMX*NATOMY][3],dec, ydist1, ydist2, ydist, pb = 0.0;
int i,j,natom,x,y,ntube,n=10,m=10,p,q,jmax,ntube1,ntube2,num1,num2,num3,count,k,l,kc;
char ch1,ch2,ch3;
int nb[]={1, 2, 3, 5, 8, 14, 15, 17, 20, 42, 43, 45, 48, 86, 87, 89, 92, 141, 143, 146};
int im[3000][4], cnt, an[30000], dh[24000][4], dhtmp[4],bnd[3000][2], ancnt;


potcut=4*(-pow(cut,-6)+pow(cut,-12));

natom=0;

/* generating graphite structure */

for(j=0;j<NATOMY/4;j++)
{ for(i=0;i<NATOMX;i++)
  {  
    atoms[natom][0]=i*(sqrt(3.))*a;
    atoms[natom][1]=j*3.*a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a;
    atoms[natom][1]=j*3.*a+a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a+a*(sqrt(3.))/2.;
    atoms[natom][1]=j*3.*a+1.5*a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a+a*(sqrt(3.))/2.;
    atoms[natom][1]=j*3.*a+2.5*a;
    natom++;
  }   
}


/* 
for(i=0;i<natom;i++)
  printf("%f\t%f\n",atoms[i][0],atoms[i][1]);
*/


/* rotation of the structure */

teta=atan(sqrt(3.)*m/(m+2.*n));  /* chiral angle */

/* rotation point */
rot[0][0]=atoms[0][0];
rot[0][1]=atoms[0][1];
for(i=1;i<natom;i++)
{
  x=atoms[i][0];
  y=atoms[i][1];

  phi=sqrt(x*x+y*y);
  alpha=acos(x/phi);

  rot[i][0]=phi*cos(alpha+teta);
  rot[i][1]=phi*sin(alpha+teta);

}


/*
for(i=0;i<natom;i++)
  printf("%f\t%f\n",rot[i][0],rot[i][1]);
*/

/* selection of atoms for the nanotube */

cir=sqrt(3.)*a*sqrt(m*m+m*n+n*n); /* length of the chiral vector */
xmin=10.1*(sqrt(3.))*a;
xmax=xmin+cir;
ymin=39.1*3.*a;
ymax=43.1*3.*a;

ntube=0;

for(i=0;i<natom;i++)
{ 
  x=rot[i][0];
  y=rot[i][1];
  if((x<xmax) && (x>xmin) && (y<ymax) && (y>ymin))
  { sheet[ntube][0]=x;
    sheet[ntube][1]=y;
    ntube++;
  }
}

/*
for(i=0;i<ntube;i++)
  printf("%f\t%f\n",sheet[i][0],sheet[i][1]);
*/



/* projection into a tube after changing the origin of the axis  */

r=cir/(2.*M_PI);
for(i=0;i<ntube;i++)
{ tube1[i][0]=r*cos(sheet[i][0]/r-xmin);
  tube1[i][1]=sheet[i][1]-(ymax+ymin)/2;
  tube1[i][2]=r*sin(sheet[i][0]/r-xmin);
}

ntube1=ntube;




n=15;
m=15;

natom=0;

/* generating graphite structure */

for(j=0;j<NATOMY/4;j++)
{ for(i=0;i<NATOMX;i++)
  {  
    atoms[natom][0]=i*(sqrt(3.))*a-a*sqrt(3.)/2.;
    atoms[natom][1]=j*3.*a-0.5*a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a-a*sqrt(3.)/2;
    atoms[natom][1]=j*3.*a+a-0.5*a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a+a*(sqrt(3.))/2.-a*sqrt(3.)/2;
    atoms[natom][1]=j*3.*a+1.5*a-0.5*a;
    natom++;
    atoms[natom][0]=i*(sqrt(3.))*a+a*(sqrt(3.))/2.-a*sqrt(3.)/2;
    atoms[natom][1]=j*3.*a+2.5*a-0.5*a;
    natom++;
  }   
}


/* 
for(i=0;i<natom;i++)
  printf("%f\t%f\n",atoms[i][0],atoms[i][1]);
*/


/* rotation of the structure */

teta=atan(sqrt(3.)*m/(m+2.*n));  /* chiral angle */

/* rotation point */
rot[0][0]=atoms[0][0];
rot[0][1]=atoms[0][1];

for(i=1;i<natom;i++)
{
  x=atoms[i][0];
  y=atoms[i][1];

  phi=sqrt(x*x+y*y);
  alpha=acos(x/phi);

  rot[i][0]=phi*cos(alpha+teta);
  rot[i][1]=phi*sin(alpha+teta);

}


/*
for(i=0;i<natom;i++)
  printf("%f\t%f\n",rot[i][0],rot[i][1]);
*/

/* selection of atoms for the nanotube */

cir=sqrt(3.)*a*sqrt(m*m+m*n+n*n); /* length of the chiral vector */
xmin=10.1*(sqrt(3.))*a;
xmax=xmin+cir;
ymin=36.1*3.*a;
ymax=46.1*3.*a;


ntube=0;

for(i=0;i<natom;i++)
{ 
  x=rot[i][0];
  y=rot[i][1];
  if((x<xmax) && (x>xmin) && (y<ymax) && (y>ymin))
  { sheet[ntube][0]=x;
    sheet[ntube][1]=y;
    ntube++;
  }
}

/*
for(i=0;i<ntube;i++)
  printf("%f\t%f\n",sheet[i][0],sheet[i][1]);
*/



/* projection into a tube after changing the origin of the axis */

r=cir/(2.*M_PI);
for(i=0;i<ntube;i++)
{ tube2[i][0]=r*cos(sheet[i][0]/r-xmin);
  tube2[i][1]=sheet[i][1]-(ymax+ymin)/2.;
  tube2[i][2]=r*sin(sheet[i][0]/r-xmin);
}

ntube2=ntube;


dec=-tube2[0][1];


for(i=0;i<ntube1;i++)
{ 
  tube1[i][0]+=r;
  tube1[i][1]+=dec;
  tube1[i][2]+=r;
}

for(i=0;i<ntube1;i++)
{ 
  tube1[i][0]=tube1[i][0]/100;
  tube1[i][1]=tube1[i][1]/100;
  tube1[i][2]=tube1[i][2]/100;
}



for(i=0;i<ntube2;i++)
{ 
  tube2[i][0]+=r;
  tube2[i][1]+=dec;
  tube2[i][2]+=r;
}

for(i=0;i<ntube2;i++)
{ 
  tube2[i][0]=tube2[i][0]/100;
  tube2[i][1]=tube2[i][1]/100;
  tube2[i][2]=tube2[i][2]/100;
}


// print  all atoms 
/*
for(i=0;i<ntube1;i++)
  printf("%f\t%f\t%f\n",tube1[i][0],tube1[i][1],tube1[i][2]);


for(i=0;i<ntube2;i++)
  printf("%f\t%f\t%f\n",tube2[i][0],tube2[i][1],tube2[i][2]);
*/



for(i=0;i<ntube2;i++)
{  
  tube1[i+ntube1][0] = tube2[i][0];  
  tube1[i+ntube1][1] = tube2[i][1];
  tube1[i+ntube1][2] = tube2[i][2];
}
ntube1+=ntube2;
   
if(argc>1)
    if (argv[1][0] == '-') 
	if ((argv[1][1] == 'p'))
{
    

    // atom pdb tube 1



    for(i=0;i<ntube1;i++)
	{ 
	    num1=(i+1)/676;
	    num2=(i+1-676*num1)/26;
	    num3=((i+1) % 26);
	    
	    ch1='C'+num1;
	    ch2='A'+num2;
	    ch3='A'+num3;
	    
	    if (i < ntube1-ntube2)
		printf("ATOM  %5d  %c%c%c NTI     1     %7.2f %7.2f %7.2f  1.00 20.00      1CRN\n",i+1,ch1,ch2,ch3,tube1[i][0],tube1[i][1],tube1[i][2]);
	    else
		printf("ATOM  %5d  %c%c%c NTO     1     %7.2f %7.2f %7.2f  1.00 20.00      1CRN\n",i+1,ch1,ch2,ch3,tube1[i][0],tube1[i][1],tube1[i][2]);

	} 
    printf("END\n");
    exit(0);
} else if (argv[1][1]=='c') 
    pb = atof(argv[2]);


//ntube1 = i;




    // atom psf tube 1

printf("PSF\n\n%8d !NTITLE\n REMARKS FILENAME=\"cnt.psf\"\n",3);
printf(" REMARKS Autogenerated carbon nanotube structure file\n");
printf(" REMARKS written by Mark Dequesne and Morris Chukhman\n\n");

printf("%8d !NATOM\n",ntube1);
for(i=0;i<ntube1;i++)
{ 
    num1=(i+1)/676;
    num2=(i+1-676*num1)/26;
    num3=((i+1) % 26);
 
    ch1='C'+num1;
    ch2='A'+num2;
    ch3='A'+num3;

    if (i < ntube1-ntube2)
	printf("   %5d 1CRN 1    NTI  %c%c%c  CX1    0.000000E+00   12.0110           0\n",i+1,ch1,ch2,ch3);
    else
	printf("   %5d 1CRN 1    NTO  %c%c%c  CX1    0.000000E+00   12.0110           0\n",i+1,ch1,ch2,ch3);

    
}


//ntube1+=ntube2;




// conect psf tube 1

/* initialize angle array */
  for(i=0; i<3000;i++)
      for(j=0;j<4;j++)
	  im[i][j]=0;

  count=0;
  for(i=0;i<ntube1;i++)
      for(j=i+1;j<ntube1;j++)
{
    
    ydist1=abs(tube1[i][1]-tube1[j][1]);
    ydist2=abs(abs(tube1[i][1]-tube1[j][1])-0.0);
    ydist = ( ydist1 < ydist2 ) ? ydist1 : ydist2;
    rad=sqrt(pow(tube1[i][0]-tube1[j][0],2)+pow(ydist,2)
	     +pow(tube1[i][2]-tube1[j][2],2));
    
    kc=0;
    
    //for(k=0;k<20;k++)
    //if(nb[k]==(i+1)) kc++;

    //for(k=0;k<20;k++)
    //if(nb[k]==(j+1)) kc++;

    //    for(k=0;(k<20)&&((i+1)!=nb[k]);k++);
	//   printf("%d\n",nb[10]);
    //if (k==20)         //((i+1)!=n[k])
  //  for(l=0;(l<20)&&((j+1)!=nb[l]);l++);

    
   //    if ( ((l!=20) || (k!=20) ))           //   ((j+1)!=n[l])

    
            if ((rad<1.6) && (kc<2))
            { count++;
	    //       printf("%8d%8d",i+1,j+1);
	    //if ((count % 4)==0) printf("\n");

		/* update improper array */
	    im[i+1][++im[i+1][0]] = j+1;
	    im[j+1][++im[j+1][0]] = i+1;
	//	printf("%d %d %d\n",na[i+1][0], i+1, j+1);
		bnd[count][0] = i+1;
		bnd[count][1] = j+1;

            }
}

for (cnt=0;cnt<count;cnt++)
{
    j = 0;
    for (i=1; i<4; i++)
	if(im[bnd[cnt][0]][i] != bnd[cnt][1])
	    dhtmp[j++] = im[bnd[cnt][0]][i];
    
    for (i=1; i<4; i++)
	if(im[bnd[cnt][1]][i] != bnd[cnt][0])
	   dhtmp[j++] = im[bnd[cnt][1]][i];
	   
    for(i=0; i < 2; i++)
	for(j=0; j<2; j++)
	    //	    for(k=0; k<4; k++)
	    {
		dh[cnt*4+i*2+j][0] = dhtmp[i];
		dh[cnt*4+i*2+j][1] = bnd[cnt][0];
		dh[cnt*4+i*2+j][2] = bnd[cnt][1];
		dh[cnt*4+i*2+j][3] = dhtmp[2+j];
	    }
}
	 

/* fill up angle array */
for (cnt=0; cnt<ntube1; cnt++)
{ 
    //   for (i = 0; ( (i < 9) && im[cnt+i][0]; i++);
    if (im[cnt][0] > 1) { 
	an[ancnt++] = im[cnt][1]; 
	an[ancnt++] = cnt;
	an[ancnt++] = im[cnt][2];
    }
    if (im[cnt][0] > 2) {
	an[ancnt++] =  im[cnt][1]; 
	an[ancnt++] = cnt;
	an[ancnt++] = im[cnt][3]; 
	an[ancnt++] = im[cnt][2]; 
	an[ancnt++] = cnt;
	an[ancnt++] = im[cnt][3];
	//ancnt++;
    }
}

printf("\n%8d !NBOND: bonds\n",count);
for (cnt=1;cnt<count;cnt+=4) {
    for(i=0;i<4;i++) 
	printf("%8d%8d",bnd[cnt+i][0],bnd[cnt+i][1]);
    printf("\n");
}


printf("\n%8d !NTHETA: angles\n",ancnt/3);
for (cnt=0; cnt<ancnt; cnt++)
{ 
              printf("%8d",an[cnt]);
              if ((cnt % 9)==8) printf("\n");
}

count = 1;
printf("\n%8d !NPHI: dihedrals\n",count*4-4);
for(cnt=4; cnt<count*4; cnt+=2) { 
    printf("%8d%8d%8d%8d",dh[cnt][0],dh[cnt][1],dh[cnt][2],dh[cnt][3]);
    printf("%8d%8d%8d%8d\n",dh[cnt+1][0],dh[cnt+1][1],dh[cnt+1][2],dh[cnt+1][3]);
}
    
ntube1 = 0;
printf("\n%8d !NIMPHI: impropers\n",ntube1);
for(cnt=1; cnt < ntube1; cnt++) { 
    printf("%8d%8d%8d%8d",cnt,im[cnt][1],im[cnt][2],im[cnt][3]);
    printf("%8d%8d%8d%8d\n",++cnt,im[cnt][1],im[cnt][2],im[cnt][3]);
}

printf("\n%8d !NDON: donors\n",0);

printf("\n%8d !NACC: acceptors\n",0);

printf("\n%8d !NNB\n",0);

printf("\n%8d%8d !NGRP\n",1,0);

}