// Code to generate lamps input file including atom sorting into three regions. 
// This code converts x y z coordinates to lammps input format
// and assign upper boundary and lower boundary atoms into atom number 1~N1, and// atom number 1~N2

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

int main(int argc,char *argv[])
{
 // void skipInitial(FILE *tensor,int n);
  FILE *fp1, *fp2, *fp3;
  int  i,j,k;
  int Natom,Nframe;  
  int Nfreq;
  double x,y,z,xlo,xhi,ylo,yhi,zlo,zhi;
  double **CreateMatrix(int rows,int cols);

  printf("usage: ./a.out in.dat out.dat atom# box_x box_y box_z z1 z2 (boundary coordinates)\n");
  if(argc!=9)
  {
    printf("Check the argument of your command\n");
    exit(0);
  }

  fp1   = fopen(argv[1],"r");
  fp2   = fopen(argv[2],"w");
  Natom =atoi(argv[3]);
  xhi   = atof(argv[4]);
  yhi   = atof(argv[5]);
  zhi   = atof(argv[6]);
  z1    = atof(argv[7]);
  z2    = atof(argv[8]);

  int    N_act,N_lower,N_upper;
  double **act, **lower, **upper;

  act = CreateMatrix(Natom,3);
  lower = CreateMatrix(Natom,3);
  upper = CreateMatrix(Natom,3);
  
  Nact=0;
  Nlower=0;
  Nupper=0;
  
  fprintf(fp2,"LAMMPS\n\n");
  fprintf(fp2,"     %d  atoms\n\n",Natom);
  fprintf(fp2,"     %d  atom types \n\n",1);
  fprintf(fp2,"     %5.3f %5.3f xlo xhi\n", 0.0, xhi);
  fprintf(fp2,"     %5.3f %5.3f ylo yhi\n", 0.0, yhi);
  fprintf(fp2,"     %5.3f %5.3f zlo zhi\n\n", 0.0, zhi);
  fprintf(fp2,"     Masses\n\n");
  fprintf(fp2,"     %d  %f\n\n", 1, 58.5934);//for Nickel
  fprintf(fp2,"     Atoms\n\n"); 

  for(k=0;k<Natom;k++)
  { 
    fscanf(fp1,"%lf %lf %lf \n",&x,&y,&z);

    if (-1 <=  z  <=  z1){
         
         lower[Nlower][0]=x;
         lower[Nlower][1]=y;
         lower[Nlower][2]=z;
         Nlower++;
    }
    else if ( z2 <= z <= zhi+1){
         upper[Nupper][0]=x;
         upper[Nupper][1]=y;
         upper[Nupper][2]=z;
         Nupper++;
    }  
    else{
         act[Nact][0]=x;
         act[Nact][1]=y;
         act[Nact][2]=z;
         Nact++;
    }
        
  }

  if(Nlower + Nact + Nupper != Natom){
     printf("error in sorting boundary atoms! Nlower:%d Nupper:%d Nact:%d\n",Nlower,Nact,Nupper);
     exit (0);
  }

  for (i=0;i<Nlower;i++)
     fprintf(fp2," %d  %d %f %f %f \n", i+1, 1, lower[Nlower][0],lower[Nlower][1],lower[Nlower][2]);
 
  for (i=0;i<Nupper;i++)
     fprintf(fp2," %d  %d %f %f %f \n", i+Nlower+1, 1, upper[Nupper][0],upper[Nupper][1],upper[Nupper][2]);
  
  for (i=0;i<Nact;i++)
     fprintf(fp2," %d  %d %f %f %f \n", i+Nlower+Nupper+1, 1, upper[Nupper][0],upper[Nupper][1],upper[Nupper][2]);

     printf("lower boundary atoms:1~%d, upper boundary atoms:%d~%d\n ", Nlower,Nlower+1,Nupper+Nlower);                                                                                                 

  
 
  fclose(fp1);
  fclose(fp2);

 
}

double ** CreateMatrix(int rows,int cols)
{
   int  i;
   double  **m;
                                                                                
   m = calloc((unsigned int) rows,sizeof(double *));
   for (i=0; i < rows; i++) {
      m[i] = calloc((unsigned int) cols,sizeof(double ));
   }
                                                                                
   return m;
}