#include<iostream>
#include<fstream>
#include<vector>
#include<cstdlib>
using namespace std;

int main (int argc, const char *argv[]) {
   int Nch = atoi(argv[1]);
   int N = atoi(argv[2]);
   int Nframes = atoi(argv[3]);
   double halfbox = atof(argv[4]); // half-box size; currently works only for cubic cells
   const char *inpos = argv[5]; // Name of the LAMMPS config-dump file
   const char *inbon = argv[6]; // Name of the LAMMPS bond-dump file

   ifstream fin1, fin2;
   fin1.open(inpos);  fin2.open(inbon);
   char ofname[50];
   ofstream fout;
   sprintf(ofname,"configs.Nch%d.N%d.reorg",Nch,N);
   fout.open(ofname);
   ofstream fout2;
   fout2.open("nextatomlist");
   ofstream fout3;
   sprintf(ofname,"data.Nch%d.N%d.frame.%d.reorg",Nch,N,Nframes);
   fout3.open(ofname);

   int f, i, j;
   int dum, type, nextAtomID, lastAtomID;
   char sponge1[255];   char sponge2[255];  
   int index1, index2;
   vector <int> nextatom;
   nextatom.resize(Nch*N);
   for (i = 0; i < Nch*N; i++) nextatom[i] = -1;

   vector <double> x, y, z;
   x.resize(Nch*N);  y.resize(Nch*N);  z.resize(Nch*N);

// for resetting images
   int ix, iy, iz, tid, lid;
   double DX, DY, DZ, xlo, xhi, ylo, yhi, zlo, zhi;
   double bx, by, bz;

 for (f = 0; f < Nframes; f++) {
    cout << "Frame " << f << endl;

// get input file headers
    for (i = 0; i < 9; i++) fin2.getline(sponge2,255);
    for (i = 0; i < 5; i++) { 
        fin1.getline(sponge1,255);
        fout << sponge1 << endl;
    }
    fin1 >> xlo >> xhi >> ylo >> yhi >> zlo >> zhi;
    DX = xhi - xlo;	DY = yhi - ylo;		DZ = zhi - zlo;
    fout << xlo << " " << xhi << endl << ylo << " " << yhi << endl << zlo << " " << zhi << endl;
    for (i = 0; i < 2; i++) {
       fin1.getline(sponge1,255);
       if (i) fout << sponge1 << endl; 
    }

// get bonds from end of DBH run
   for (i = 0; i < Nch*(N-1); i++) {
     fin2 >> dum >> dum >> index1 >> index2;
     index1 -= 1;  index2 -= 1;
     nextatom[index1] = index2;
   }
   fin2.getline(sponge2,255);

// get coords: reset images
   for (i = 0; i < Nch*N; i++) {
      fin1 >> index1 >> dum >> dum;
      index1 -= 1;  
      fin1 >> x[index1] >> y[index1] >> z[index1] >> dum >> dum >> dum;
   }
   fin1.getline(sponge1,255);

// special for data file from last frame
 if (f == Nframes-1) {
   fout3 << "LAMMPS FENE chain data file" << endl << endl;
   fout3 << Nch*N << " atoms" << endl << Nch*(N-1) << " bonds" << endl << Nch*(N-2) << " angles" << endl;
   fout3 << "2 atom types" << endl << "1 bond types" << endl << "1 angle types" << endl << endl;
   fout3 << -halfbox << " " << halfbox << " xlo xhi" << endl;
   fout3 << -halfbox << " " << halfbox << " ylo yhi" << endl;
   fout3 << -halfbox << " " << halfbox << " zlo zhi" << endl;
   fout3 << endl << "Masses" << endl << endl << "1 1.0" << endl << "2 1.0" <<  endl << endl << "Atoms" << endl << endl;
   }

// Reordering into molecules according to new bond list
// make chain-end atoms type 2 for PPA
   for (i = 0; i < Nch; i++) {
      ix = 0; iy = 0; iz = 0;
      for (j = 0; j < N; j++) {
         if (j == 0 || j == N-1) type = 2;   else type = 1;
         fout << i*N + j + 1 << " " << i+1 << " " << type << " ";
         if (f == Nframes-1) fout3 << i*N + j + 1 << " " << i+1 << " " << type << " ";
         if (j == 0) nextAtomID = i*N + j;
         else { lastAtomID = nextAtomID; nextAtomID = nextatom[nextAtomID]; }

// reset images
         if (j > 0) {
            tid = nextAtomID;
            lid = lastAtomID;
            bx = x[tid] - x[lid];
            by = y[tid] - y[lid];
            bz = z[tid] - z[lid];
            if (bx < -DX/2) ix++;
            else if (bx > DX/2) ix--;
            if (by < -DY/2) iy++;
            else if (by > DY/2) iy--;
            if (bz < -DZ/2) iz++;
            else if (bz > DZ/2) iz--;
         }
         
         fout << x[nextAtomID] << " " << y[nextAtomID] << " " << z[nextAtomID] << " " << ix << " " << iy << " " << iz << endl;
         if (f == Nframes-1) fout3 << x[nextAtomID] << " " << y[nextAtomID] << " " << z[nextAtomID] << " " << ix << " " << iy << " " << iz << endl;
      }
      for (j = 0; j < N; j++) fout2 << i*N + j + 1 << " " << nextatom[i*N+j] << endl;
    }
   
// bonds
if (f == Nframes-1) {
  fout3 << endl << "Bonds" << endl << endl;
  for (i = 0; i < Nch; i++) {
      for (j = 0; j < N-1; j++) {
         fout3 << i*(N-1) + j + 1 << " 1 " << i*N + j + 1 << " " << i*N + j + 2 << endl;
      }
  }

  fout3 << endl << "Angles" << endl << endl;
  for (i = 0; i < Nch; i++) {
     for (j = 0; j < N-2; j++) {
        fout3 << i*(N-2) + j + 1 << " 1 " << i*N + j + 1 << " " << i*N + j + 2 << " " << i*N + j + 3 << endl;
     }
    }
  }
}
   fin1.close();  fin2.close();
   fout.close();  fout2.close();
 
   return 0;

}