using namespace std;
#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<fstream>
#include<vector>

int main(){
double Lx,Ly,Lz;
int nbead,nfiber,ntype,nbond,nangle,type,mol;
int nbeadpoly,ntypepol;
cout <<"How many poly: "<<endl;
cin>> nfiber;
cout << "How many atom types: "<<endl;
cin >> ntypepol;
cout << "How many bead per poly: "<<endl;
cin >> nbeadpoly;
cout << "Box size Lx \n Ly \n Lz: "<<endl;
//cin >> Lx;
//cin>>Ly;
//cin >> Lz;
int closed;
cout << "Open(0) or Closed(1)?" <<endl;
cin >> closed;
Lx=Ly=Lz=40;
int nmax = int(nfiber*nbeadpoly + 0.5*nbeadpoly); 

vector<int> atom_type(nmax);
vector<int> mol_type(nmax);
vector<int> nangle1_type(nmax);
vector<int> bond_type(nmax);
vector<int> Nfib(nfiber);
vector<vector< int> > chainbead(nbeadpoly,vector<int>(nfiber));
//int chain[Nmax];
double x,y,z;
vector<vector<int> > angle(nmax,vector<int>(4));
vector<vector<int> > angle1(nmax,vector<int>(3));
vector<vector<int> > bond(nmax,vector<int>(2));
vector<vector<double> >PosFib(nfiber,vector<double>(3));
vector<vector<vector< double> > >position(nbeadpoly, PosFib);

nbead =0;
ofstream write;
ifstream read;
read.open("initial_configuration.dat");
write.open("trefoil_lammps_patchy.out");
//int count=0;
for(int j=0;j<nfiber;j++){ 
	for(int i=0;i<nbeadpoly;i++){
	read >> mol >> type>> position[i][j][0] >> position[i][j][1] >>position[i][j][2] ;
	atom_type[nbead]=type;
	mol_type[nbead]=mol;
	chainbead[i][j]=nbead;
//	polymer[i][j][0] = x;
//	polymer[i][j][1] = y;
//	polymer[i][j][2] = z;	
//	cout << nbead << " " <<type <<" "<< polymer[i][j][0] << " "<< polymer[i][j][1]<<" "<<polymer[i][j][2]<<endl;
nbead++;
	}
}


//while(read.eof()){
//nbead++;
//read >> type >> x >> y >> z;
//cout << nbead <<" "<< x<<endl; 
//}


read.close();

nbond=0;
nangle=0;
int nangle1=0;

int nbeadT=0;
//EVEN BEADS (counting from 0)
int endBonds;
if(closed==1) endBonds=nbeadpoly-1;
else endBonds=nbeadpoly-2;
for(int j=0;j<nfiber;j++){
	for(int i =0;i<endBonds;i++){
	nbeadT++;
	if(i%2==0){
		if((nbeadT+1)%nbeadpoly==0) {		//LAST BOND: M-2 ---- 0 and M-2 ---- M-1 !!!!!!!
		//bond[nbond][0] = chainbead[i][j];
		//bond[nbond][1] = chainbead[i+1][j];
		//bond_type[nbond] = 2;
		//nbond=nbond+1;
		bond[nbond][0] = chainbead[i][j];		
		bond[nbond][1] = chainbead[0][j];
		bond_type[nbond]=1;
		}
		if((nbeadT+1)%nbeadpoly!=0){
		//bond[nbond][0] = chainbead[i][j];
		//bond[nbond][1] = chainbead[i+1][j]; NO SIDEWAY BOND IT's RIGID BODY!!
		//bond_type[nbond] = 2;
		//nbond=nbond+1;
		bond[nbond][0] = chainbead[i][j];		
		bond[nbond][1] = chainbead[i+2][j];
		bond_type[nbond] = 1;
		}
	nbond++;
	}
	}
}
//ODD BEADS(counting from 0)
//for(int j=0;j<nfiber;j++){
//	for(int i =0;i<nbeadpoly;i++){ //Qui e`  -1 o 0 
//	nbeadT++;
//	if(i%2==1){
//		if((nbeadT)%nbeadpoly==0) {		//LAST M-1 --- 1 !!!!!!!
//		bond[nbond][0] = chainbead[i][j];
//		bond[nbond][1] = chainbead[1][j];
//		bond_type[nbond] = 1;
//		}
//		if((nbeadT)%nbeadpoly!=0){
//		bond[nbond][0] = chainbead[i][j];
//		bond[nbond][1] = chainbead[i+2][j];
//		bond_type[nbond] = 1;
//		}
//	nbond++;
//	}
//	}
//}	
cout << "nbonds " << nbond<<endl;
//ANGLES
int endAngle;
if(closed==1) endAngle=nbeadpoly-3;
else endAngle=nbeadpoly-4;
for(int j=0;j<nfiber;j++){ 
	for(int i =0;i<endAngle;i++){ //qui e` -3 o -4  nbeadpoly-3
	if(i%2==0){		//EVEN
	if(i==nbeadpoly-4){		//LAST ANGKE M-4 --- M-2 --- 0 
	angle1[nangle1][0] = chainbead[i][j];
	angle1[nangle1][1] = chainbead[i+2][j];
	angle1[nangle1][2] = chainbead[0][j];
	}
	if(i!=nbeadpoly-4){
	angle1[nangle1][0] = chainbead[i][j];
	angle1[nangle1][1] = chainbead[i+2][j];
	angle1[nangle1][2] = chainbead[i+4][j];
	}
	nangle1_type[nangle1]=1;
	nangle1=nangle1+1;
	}
	}
}
//ANGLES 90 degress Patch2-Bead1-Bead3
int endAngle2;
if(closed==1) endAngle2=nbeadpoly;
else endAngle2=nbeadpoly-1;
for(int j=0;j<nfiber;j++){ 
	for(int i =0;i<endAngle2;i++){ //qui e` -2 o -3
	if(i%2==1){		//ODD
	if(i==nbeadpoly-1){		//LAST ANGLE M-3 --- M-1 --- 1 
	angle1[nangle1][0] = chainbead[i][j];
	angle1[nangle1][1] = chainbead[i-1][j];
	angle1[nangle1][2] = chainbead[0][j];
	}
	if(i!=nbeadpoly-1){
	angle1[nangle1][0] = chainbead[i][j];
	angle1[nangle1][1] = chainbead[i-1][j];
	angle1[nangle1][2] = chainbead[i+1][j];
	}
	nangle1_type[nangle1]=2;
	nangle1=nangle1+1;
	}
	}
}

cout << "nangles1 " << nangle1<<endl;

//DIHEDRAL 
int endDihe;
if(closed==1) endDihe=nbeadpoly-1;
else endDihe=nbeadpoly-2;
for(int j=0;j<nfiber;j++){
	for(int i =0;i<endDihe;i++){ //qui e` -1 o -2
	if(i%2==0){
	if(i==nbeadpoly-2){		//LAST DIHEDRAL M-1 --- M-2 --- 0 --- 1
	angle[nangle][0] = chainbead[i+1][j];
	angle[nangle][1] = chainbead[i][j];
	angle[nangle][2] = chainbead[0][j];
	angle[nangle][3] = chainbead[1][j];	
	}
	if(i!=nbeadpoly-2){
	angle[nangle][0] = chainbead[i+1][j];
	angle[nangle][1] = chainbead[i][j];
	angle[nangle][2] = chainbead[i+2][j];
	angle[nangle][3] = chainbead[i+3][j];	
	}
	nangle=nangle+1;
	}
	}
}
cout << "nangle " << nangle<<endl;

write << " LAMMPS data file from restart file: timestep = 0,\t procs = 1"<<endl;
write << nbead << " atoms "<<endl;
write <<  nbond << " bonds "<<endl;
write <<  nangle1 << " angles "<<endl;
write << nangle << " dihedrals "<<endl;
write << "\n";
write << ntypepol << " atom types "<<endl;
write << 1 << " bond types "<<endl;
write << 2 << " angle types "<<endl;
write << 1 << " dihedral types "<<endl;
write << "\n";
write << -Lx/2.0 << " " << Lx-Lx/2.0 << " xlo xhi"<<endl;
write << -Ly/2.0 << " " << Ly-Ly/2.0 << " ylo yhi"<<endl;
write << -Lz/2.0 << " " << Lz-Lz/2.0 << " zlo zhi"<<endl;
write << "\n Masses \n"<<endl;

for(int j=1; j<=ntypepol;j++) write << j << " " << 1 << endl; 
 
write << "\n Atoms \n"<<endl;

int nbeadTemp=0;
for(int j=0;j<nfiber;j++){
	for(int i=0;i<nbeadpoly;i++){
	write << nbeadTemp+1<<" "<<mol_type[nbeadTemp]<<" "<<atom_type[nbeadTemp]<<"\t "<<position[i][j][0] << " " << position[i][j][1] << " " << position[i][j][2] << " \t" << 0 << " " << 0<< " " << 0<<endl;
nbeadTemp++;
	}
}

write << "\n Velocities \n" <<endl;
 
for(int j=0;j<nbead;j++) write<<j+1 << " "<<0 << " "<<0<<" "<<0 <<endl;

write << "\n Bonds \n"<<endl; 
for(int i=0;i<nbond;i++) {
write << i+1 <<" "<< bond_type[i] <<" "<< bond[i][0]+1<<" "<<bond[i][1]+1 << endl;
}

write << "\n Angles \n"<<endl; 
for(int i=0;i<nangle1;i++) write << i+1 <<" "<< nangle1_type[i] <<" "<< angle1[i][0]+1<<" "<<angle1[i][1]+1 <<" "<< angle1[i][2]+1 << endl;


write << "\n Dihedrals \n"<<endl; 
for(int i=0;i<nangle;i++) write << i+1 <<" "<< 1 <<" "<< angle[i][0]+1<<" "<<angle[i][1]+1 <<" "<< angle[i][2]+1<< " " << angle[i][3]+1 << endl;

write.close();
return 0;
}