// This code analyzing shear stress and shear strain from the ouput file of lammps containg stress tensor per atom

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

int main(int argc,char *argv[])
{
  void skipInitial(FILE *tensor,int n);
  FILE *fp1, *fp2;
  int  i,j;
  double f, stress,strain,plate_v,height,stepSize;
  int Natom,Nframe;  
  int Nfreq;
  double sxx,syy,szz,sxy,sxz,syz;

  printf("usage: ./a.out in.dat out.dat atom# frame# freq velocity(upper plate)(A/ps) timestep box_height\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]);
  Nframe = atoi(argv[4]);
  Nfreq = atoi(argv[5]);
  plate_v = atof(argv[6]);
  stepSize = atof(argv[7]);
  height   = atof(argv[8]);
 
 for (j=0;j<Nframe;j++){
  skipInitial(fp1,9);
  f=0;
  for(i=0;i<Natom;i++)
  {
    fscanf(fp1,"%lf %lf %lf %lf %lf %lf\n",&sxx,&syy,&szz,&sxy,&sxz,&syz);
    f=f + sxz;
  }
  stress=f/Natom;
  strain=(plate_v*Nfreq*stepSize*j)/height;
 // printf("%e %d %e %d %e\n",plate_v,Nfreq,stepSize,j,height);
  fprintf(fp2,"%e %e \n",strain,stress);
  }
  fclose(fp1);
  fclose(fp2);

 
}

void skipInitial(FILE *tensor,int n)
{
  int i;
  for(i=0;i<n;i++)
    {
      do
	{
	}
      while(fgetc(tensor)!='\n');
    }

  return;
}