//######################################################################################//
//This version incorporates a randomisation routine for atom allocation overflow		//
//######################################################################################//

#define VERSION 080222
#include <stdio.h>
#include <math.h>
#include <stddef.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

//---qq for debuging---//
#define qq puts("bla\n");

#define NR_END 1
#define FREE_ARG char*

void nrerror(char[]);
float *vector(long,long);
int *dvector(int,int);
void free_vector(float*,long,long);
void Randomize(void);
int Random(int);


int main(int argc,char *argv[])
{
//---Variable declaration---//
int i,j,k,temp,natoms,cmp,centsetcount,replace,low,high,overflow;
int centtype,bondatom,coord,repsetcount,ntypes,test,newat,*ovrflw;
float temp1,temp2,temp3,temp4,tetra;
float rmin,rmax,posi,x,y,z;
float dx,dy,dz;
char line[100],junk,labeltype[100],yn;

//---File pointer for input data---//
FILE *in;

if(argc!=2)
        {
        fprintf(stderr,"\nYou need to supply ./subatom filename\n");
        fprintf(stderr,"The input file is a .chem3d file.\n");
        fprintf(stderr,"Output is directed to stdout, messages to stderr.\n");		
        return 1;
        }

//---Open file or dump error msg---//
        if (!(in=fopen(argv[1],"r")))
                {
                fprintf(stderr,"Unable to open file %s\n",argv[1]);
                return 1;
                }
//---first line of .chem3d file should be number of atoms in file---//
fgets(line,100,in);
sscanf(line,"%d",&temp);
natoms=temp;
//---Setup matrix for input data---//
float data[natoms+1][6];

//---Reads data in from file---//
i=1;j=0;
cmp=junk;
while(i!=natoms+1){
	fgets(line,100,in);
	sscanf(line,"%s%d%f%f%f",&junk,&temp,&dx,&dy,&dz);
	data[i][1]=i;
	data[i][3]=dx;
	data[i][4]=dy;
	data[i][5]=dz;
//---Checks number of atoms and atom type---//
	if(junk!=cmp) j++;
	cmp=junk;
	data[i][2]=j;
	i++;
	}

ntypes=j;
fprintf(stderr,"\nTHIS PROGRAM HAS NOT BEEN BUG FIXED FOR STUPID USER INPUTS!!!! - BETA USE ONLY\n");
fprintf(stderr,"\nThere are %d atoms and %d atom types.\n",natoms,ntypes);
//---Ask user for central atom of tetrahedron---//
fprintf(stderr,"Central Atoms type: ");
scanf("%d",&centtype);
//--Check to see how many atoms of type central there are in the file---//
j=0;
for(i=1;i<=natoms;i++) if (data[i][2]==centtype) j++;
fprintf(stderr,"Number of atoms of type %d: %d\n",centtype,j);
//---Move these atoms to another matrix to reduce computation time during coordination run---//
float centset[j+1][6];
centsetcount=1;
for(i=1;i<=natoms+1;i++){
	if(data[i][2]==centtype){
		centset[centsetcount][1]=data[i][1];
		centset[centsetcount][2]=data[i][3];
		centset[centsetcount][3]=data[i][4];
		centset[centsetcount][4]=data[i][5];
		centset[centsetcount][5]=0;
		centsetcount++;
	}
}
centsetcount=centsetcount-1;
//---Ask user for bond atoms type and parameters---//
fprintf(stderr,"Bonding Atom type: ");
scanf("%d",&bondatom);
fprintf(stderr,"RMin: ");
scanf("%f",&rmin);
fprintf(stderr,"RMax: ");
scanf("%f",&rmax);
fprintf(stderr,"Coordination: ");
scanf("%d",&coord);
//---Work out coordination numbers for bonding atom---//
for(i=1;i<=natoms+1;i++){
	if(data[i][2]==bondatom){
		for(j=1;j<=centsetcount+1;j++){
			x=data[i][3]-centset[j][2];
			y=data[i][4]-centset[j][3];
			z=data[i][5]-centset[j][4];
			posi=sqrt(x*x+y*y+z*z);
			if(posi>=rmin && posi<=rmax) centset[j][5]++;
		}
	}
}
//---Setting up new matrix for only those central atoms of desired coordination---//
repsetcount=0;
for(i=1;i<=centsetcount+1;i++) if(centset[i][5]==coord) repsetcount++;

float repset[repsetcount][6];
j=1;
for(i=1;i<=centsetcount+1;i++){
	if(centset[i][5]==coord){
		repset[j][1]=centset[i][1];
		repset[j][2]=centset[i][2];
		repset[j][3]=centset[i][3];
		repset[j][4]=centset[i][4];
		j++;
	}
}
//fprintf(stderr,"There are %d atoms of type %d coordinated as %d.\n",repsetcount,bondatom,coord);
//---Ask user for how many atoms to replace---//
fprintf(stderr,"Number of bonded atoms to replace: ");
scanf("%d",&replace);
//---Randomise the set, otherwise there will be clustering of replaced atoms---//
i=0;
while(i<=1000){
	low=Random(repsetcount);
	high=Random(repsetcount);
	if(low!=0 && high!=0){
		temp1=repset[low][1];
		temp2=repset[low][2];
		temp3=repset[low][3];
		temp4=repset[low][4];				
		repset[low][1]=repset[high][1];
		repset[low][2]=repset[high][2];
		repset[low][3]=repset[high][3];
		repset[low][4]=repset[high][4];						
		repset[high][1]=temp1;
		repset[high][2]=temp2;
		repset[high][3]=temp3;
		repset[high][4]=temp4;				
	}
	i++;
}
newat=ntypes+1;
i=1;
tetra=0;
test=1;
//---Replace the bonded atoms---//
for(i=1;i<=repsetcount;i++){
	for(j=1;j<=natoms;j++){
		x=repset[i][2]-data[j][3];
		y=repset[i][3]-data[j][4];
		z=repset[i][4]-data[j][5];
		posi=sqrt(x*x+y*y+z*z);
		if(posi>=rmin && posi<=rmax && test<=replace && data[j][2]==(float)bondatom){
			data[j][2]=(float)newat;
			test++;
			tetra=tetra+0.25;
			}
	}
}
//---Inform user of how many atoms are changed---//
temp1=0;
for(i=1;i<=natoms+1;i++) if (data[i][2]==(float)newat) temp1++;
fprintf(stderr,"%d atoms changed.\n",(int)temp1);
//fprintf(stderr,"%f Total/Partial Tetrahedra.\n",tetra);

//---Code to use overflow option if replace not reached.
temp2=1;
if((int)temp1!=replace){
	fprintf(stderr,"Number of atoms to change not reached. Overflow? (y/n): ");
	scanf("%s",&yn);
	if(yn=='y'){
//-----------------------Overflow code--------------------//	
		for(i=1;i<=natoms;i++){
			if(data[i][2]==bondatom) overflow++;
			}
j=1;
ovrflw=dvector(1,overflow);
		for(i=1;i<=natoms;i++){
			if(data[i][2]==bondatom){
				ovrflw[j]=data[i][1];
				j++;
				}
			}
		while(i<=1000){
		low=Random(repsetcount);
		high=Random(repsetcount);
		if(low!=0 && high!=0){
			temp1=ovrflw[low];
			ovrflw[low]=ovrflw[high];
			ovrflw[high]=temp1;
			}
		i++;
		}
		temp2=0;
		for(i=1;i<=replace-temp1;i++){
			data[ovrflw[i]][2]=(float)newat;
			temp2++;
		}
//-----------------------Overflow code--------------------//
		fprintf(stderr,"Total now changed is %d\n",(int)temp1+(int)temp2);
	}
	else{
		fprintf(stderr,"Okay but remember the amount you asked to be changed has not been reached.\n");
	}
}







//---Move replaced atoms to end of list, and make output matrix---//
float output[natoms+1][6];
j=1;
for(i=1;i<=natoms;i++){
	if(data[i][2]!=(float)newat){
		output[j][1]=data[i][1];
		output[j][2]=data[i][2];
		output[j][3]=data[i][3];
		output[j][4]=data[i][4];
		output[j][5]=data[i][5];
		j++;								
	}
}
for(i=1;i<=natoms;i++){
	if(data[i][2]==(float)newat){
		output[j][1]=data[i][1];
		output[j][2]=data[i][2];
		output[j][3]=data[i][3];
		output[j][4]=data[i][4];
		output[j][5]=data[i][5];
		j++;								
	}
}

//---Output to stdout---//
k=1;
output[0][2]=1;
fprintf(stdout,"%d\n",natoms);
fprintf(stderr,"Label for atom type 1: ");
scanf("%s",labeltype);
for(i=1;i<=natoms;i++){
	fprintf(stdout,"%s\t%d\t%f\t%f\t%f\n",labeltype,i,output[i][3],output[i][4],output[i][5]);
	if(output[i][2]!=output[i-1][2]){
		k++;
		fprintf(stderr,"Label for atom type %d: ",k);
		scanf("%s",labeltype);
		}
	}






return 0;
}

 void Randomize() {
srand( (unsigned)time( NULL ) ) ;
 }

  int Random(int Max) {
return ( rand() % Max)+ 1; 
   }

void nrerror(char error_text[])
{
    fprintf(stderr,"Numerical Recipes run-time error...\n");
    fprintf(stderr,"%s\n",error_text);
    fprintf(stderr,"...now exiting to system...\n");
    exit(1);
}

int *dvector(int nl,int nh)
{
        int *v;
        v=(int *)malloc((size_t)((nh-nl+ 1 +NR_END)*sizeof(int)));
        if (!v) nrerror("Allocation failure in dvector()");
        return v-nl+NR_END;
}
	
float *vector(long nl,long nh)
{
        float *v;
        v=(float *)malloc((size_t)((nh-nl+ 1 +NR_END)*sizeof(float)));
        if (!v) nrerror("Allocation failure in vector()");
        return v-nl+NR_END;
}

void free_vector(float *v,long nl, long nh)
{
        free((FREE_ARG)(v+nl-NR_END));
}