/*
CHANGE LOG

April 19,2007 - Initiated this program

OBJECTIVE

create a LAMMPS data file to put TIP3P water inside a simulation box

*/



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

//Random Number Functions
float RandR(int seed){
	int temp;

	srand((unsigned)seed);
	temp=rand();
	return (float)temp/(RAND_MAX);

}
//Gives Values from 0 to 1
float RandR1(void) { 
   return (rand()%RAND_MAX)/(float)RAND_MAX;
}
void RandVec3(float *p){
	float x,y,s,seed1,seed2;
	s=2.0;
	while(s>1.0){
		x=2.0 * RandR1() -1.0;
		y=2.0 * RandR1() -1.0;
		s= x * x + y * y;
	}
	p[2]=1.0 - 2.0*s;
	s=2.0 * sqrtf(1.0-s);
	p[0]=s*x;
	p[1]=s*y;
}
void RandVec3_Rectangular(float *p,float xhi, float yhi, float zhi){
                                                                                                                             
        p[0]=(2.0*RandR1()-1.0) * xhi;
        p[1]=(2.0*RandR1()-1.0) * yhi;
        p[2]=(2.0*RandR1()-1.0) * zhi;
}
//end of Random Number Functions

//Vector functions
void CrossProduct(float *A, float *B, float *C){
//AXB=C

	C[0]= A[1]*B[2] - A[2]*B[1];
	C[1]= A[2]*B[0] - A[0]*B[2];
	C[2]= A[0]*B[1] - A[1]*B[0];

}
float DotProduct(float *A, float *B){

	return (A[0]*B[0] + A[1]*B[1] + A[2]*B[2]) ;

}
float CosThetaKJ( float *A, float *B){

	float L1, L2;

	L1 = sqrtf ( A[0]*A[0] + A[1]*A[1] + A[2]*A[2] );

	L2 = sqrtf ( B[0]*B[0] + B[1]*B[1] + B[2]*B[2] );

	return (DotProduct(A,B)/ (L1*L2) );

}



void RotatePointAboutLine(float *p,float theta,float *p1,float *p2)
{

//function based on http://local.wasp.uwa.edu.au/~pbourke/geometry/rotate/example.c

	float u[3],q1[3],q2[3];
	float d,Ru;
	int k;

	/* Step 1 */
	q1[0] = p[0] - p1[0];
	q1[1] = p[1] - p1[1];
	q1[2] = p[2] - p1[2];

	u[0] = p2[0] - p1[0];
	u[1] = p2[1] - p1[1];
	u[2] = p2[2] - p1[2];

	Ru=0;
	for(k=0;k<3;k++){
		Ru+= u[k]*u[k];
	}
	Ru = sqrtf(Ru);
	
	//Normalise(&u);
	for(k=0;k<3;k++){
		u[k]/=Ru;
	}


	d = sqrtf(u[1]*u[1] + u[2]*u[2]);

	/* Step 2 */
	if (d != 0) {
		q2[0] = q1[0];
		q2[1] = q1[1] * u[2] / d - q1[2] * u[1] / d;
		q2[2] = q1[1] * u[1] / d + q1[2] * u[2] / d;
	} else {
		for(k=0;k<3;k++){
			q2[k] = q1[1];
		}
	}

	/* Step 3 */
	q1[0] = q2[0] * d - q2[2] * u[0];
	q1[1] = q2[1];
	q1[2] = q2[0] * u[0] + q2[2] * d;

	/* Step 4 */
	q2[0] = q1[0] * cosf(theta) - q1[1] * sinf(theta);
	q2[1] = q1[0] * sinf(theta) + q1[1] * cosf(theta);
	q2[2] = q1[2];

	/* Inverse of step 3 */
	q1[0] =	q2[0] * d + q2[2] * u[0];
	q1[1] =	q2[1];
	q1[2] = - q2[0] * u[0] + q2[2] * d;

	/* Inverse of step 2 */
	if (d != 0) {
		q2[0] =	q1[0];
		q2[1] =	q1[1] * u[2] / d + q1[2] * u[1] / d;
		q2[2] = - q1[1] * u[1] / d + q1[2] * u[2] / d;
	} else {
		for(k=0;k<3;k++){
			q2[k] = q1[1];
		}
	}

	/* Inverse of step 1 */
	q1[0] = q2[0] + p1[0];
	q1[1] = q2[1] + p1[1];
	q1[2] = q2[2] + p1[2];

	for(k=0;k<3;k++){
		//printf("\t%f %f %f %f\n",p[k],q1[k],p1[k],p2[k]);
		p[k]=q1[k];
	}
 
}
//End of vector functions

//gets the coordinate of oxygen and gives the coordinates of the two corresponding Hydrogens

void CreateWater(float *pOxy,float *pHyd1,float *pHyd2){

	float rOH = 0.9572;
	float theta;
	float p[3],pO[3];
	theta = 104.52 * 2.0 * M_PI / 360.0;
	//printf ("%f\n", theta);
	int k;

	for(k=0;k<3;k++){
		pO[k]=0.0;
	}

	//get coordinates of first hydrogen
	RandVec3(p);
	for(k=0;k<3;k++){
		pHyd1[k] = pO[k] + ( p[k] * rOH );
	}

	// get coordinates of second hydrogen
	RandVec3(p);
	for(k=0;k<3;k++){
		pHyd2[k] = pO[k] + ( p[k] * rOH );
	}


	float A[3],B[3],C[3];
	for(k=0;k<3;k++){
		A[k] = pHyd1[k] - pO[k];
		B[k] = pHyd2[k] - pO[k];
	}
	//get cross product
	CrossProduct(A,B,C);
	//get angle between the 2 OH Bonds
	float CheckTheta;
	CheckTheta = acosf(CosThetaKJ( A,B));
	//printf ("%f\n", CheckTheta);

	//rotate by 104.52 - CheckTheta angle
	CheckTheta = (theta - CheckTheta);
	//printf ("%f\n", CheckTheta);

	//rotate the second hydrogen to maintan an HOH angle of 104.52 deg
	RotatePointAboutLine(pHyd2,CheckTheta,pO,C);

	for(k=0;k<3;k++){
		pHyd1[k]+= pOxy[k];
		pHyd2[k]+= pOxy[k];
	}

	//Check Bondlength
	float Rab=0;
	for(k=0;k<3;k++){
		Rab+= (pHyd1[k] - pOxy[k])*(pHyd1[k] - pOxy[k]);
	}
	Rab = sqrtf(Rab);
	//printf ("%f\n", Rab);
	Rab=0;
	for(k=0;k<3;k++){
		Rab+= (pHyd2[k] - pOxy[k])*(pHyd2[k] - pOxy[k]);
	}
	Rab = sqrtf(Rab);
	//printf ("%f\n", Rab);

	//check angle
	for(k=0;k<3;k++){
		A[k] = pHyd1[k] - pOxy[k];
		B[k] = pHyd2[k] - pOxy[k];
	}

	
}

main(){
	int i,j,k;
	int nH2O;
	float pbcl[3];
	FILE *fpw1, *fpw2;
	fpw1 = fopen("Water.view","w");
	fpw2 = fopen("Water.dat","w");
	printf("\nEnter number of water molecules: ");
	scanf("%i",&nH2O);
	printf("nH2O = %i\n",nH2O);

	pbcl[0] = powf( nH2O * 29.889,1.0/3.0)/2.0;
	pbcl[1]= pbcl[0];
	pbcl[2] = pbcl [0];
	printf ("Half Boxsize is %f %f %f\n",pbcl[0],pbcl[1],pbcl[2]);

	//allocate array for water molecules
	float rH2O[nH2O][3][3];
	for(i=0;i<nH2O;i++){
		for(j=0;j<3;j++){
			for(k=0;k<3;k++){
				rH2O[i][j][k]=0;
			}
		}
	}
	
	//put water molecules in the box
	for(i=0;i<nH2O;i++){
		if(i%1000==0){
			printf("%i\n",i);
		}
		int check;
		int reps;
		reps = 0;
		while(1){
			check=0;
			float pO[3];
			RandVec3_Rectangular(pO,pbcl[0],pbcl[1],pbcl[2]);
			for(k=0;k<3;k++){
				rH2O[i][1][k]=pO[k];
			}
			//check against other Oxygen atoms
			for(j=0;j<i;j++){
				float dr[3];
				for(k=0;k<3;k++) dr[k]=0;
				float dr2;
				dr2=0;
				for(k=0;k<3;k++){
					dr[k]=(rH2O[i][1][k]-rH2O[j][1][k]);
					if(dr[k]>pbcl[k]){
						dr[k]=dr[k]-(2.0*pbcl[k]);
					}
					else if(dr[k]<-pbcl[k]){
						dr[k]=dr[k]+(2.0*pbcl[k]);
					}
					dr2+=(dr[k]*dr[k]);
				}
				dr2=sqrtf(dr2);
				if(dr2<2.0){
					check=1;
					reps++;
					break;
				}
			}
			if(check==0){
				//get the coordinates of the 2 hydrogens
				float pH1[3],pH2[3];
				CreateWater(pO,pH1,pH2);
				for(k=0;k<3;k++){
					rH2O[i][0][k] = pH1[k];
					rH2O[i][2][k] = pH2[k];
				}
				break;
			}
			else{
				if(reps==3 && i>0){
					reps=0;
					i=i-1;
					break;
				}
			}
		}
	}
	//apply periodic boundary conditions
	for(i=0;i<nH2O;i++){
		for(j=0;j<3;j++){
			for(k=0;k<3;k++){
				if(rH2O[i][j][k]>pbcl[k]){
					rH2O[i][j][k]=(rH2O[i][j][k])-((pbcl[k])*2.0);
				}
				else if (rH2O[i][j][k]< -1.0*pbcl[k] ){
					rH2O[i][j][k]=(2.0*(pbcl[k]))+ rH2O[i][j][k];
				}
			}
		}
	}
	//create viewer file
	fprintf(fpw1,"0 %f %f %f %i %i\n",pbcl[0],pbcl[1],pbcl[2],3*nH2O,2*nH2O);

	int v=0;
	for(i=0;i<nH2O;i++ ){
		for(j=0;j<2;j++){
			v++;
			fprintf(fpw1,"%i %i %i\n", v, i*(3) + j ,i*(3) + j +1 );
		}
	}
	for(i=0;i<nH2O;i++){
		fprintf(fpw1,"%f %f %f red\n",rH2O[i][0][0],rH2O[i][0][1],rH2O[i][0][2]);
		fprintf(fpw1,"%f %f %f blue\n",rH2O[i][1][0],rH2O[i][1][1],rH2O[i][1][2]);
		fprintf(fpw1,"%f %f %f red\n",rH2O[i][2][0],rH2O[i][2][1],rH2O[i][2][2]);
	}
	fclose(fpw1);

	//create LAMMPS data file
	//Build Header
	fprintf(fpw2,"LAMMPS Description\n\n");
	fprintf(fpw2,"%i atoms\n",nH2O*3);
	fprintf(fpw2,"%i bonds\n",nH2O*2);
	fprintf(fpw2,"%i angles\n\n",nH2O);
	fprintf(fpw2,"2 atom types\n");
	fprintf(fpw2,"1 bond types\n");
	fprintf(fpw2,"1 angle types\n\n");
	fprintf(fpw2,"%f %f xlo xhi\n",-1.0 * pbcl[0],pbcl[0]);
	fprintf(fpw2,"%f %f ylo yhi\n",-1.0 * pbcl[1],pbcl[1]);
	fprintf(fpw2,"%f %f zlo zhi\n\n",-1.0 * pbcl[2],pbcl[2]);
	fprintf(fpw2,"Masses\n\n");
	fprintf(fpw2,"1 1.008\n");
	fprintf(fpw2,"2 15.9994\n");

	fprintf(fpw2,"\nAtoms\n\n");
	for(i=0;i<nH2O;i++){
		for(j=0;j<3;j++){
			if(j==1){ //Oxygen
				fprintf(fpw2,"%i %i 2 -0.830 %f %f %f\n",(i*3)+j+1,i+1,rH2O[i][j][0],rH2O[i][j][1],rH2O[i][j][2]);
			}
			else{//Hydrogen
				fprintf(fpw2,"%i %i 1 0.415 %f %f %f\n",(i*3)+j+1,i+1,rH2O[i][j][0],rH2O[i][j][1],rH2O[i][j][2]);
			}
		}
	}

	fprintf(fpw2,"\nBonds\n\n");
	v=0;
	for(i=0;i<nH2O;i++ ){
		for(j=0;j<2;j++){
			v++;
			fprintf(fpw2,"%i 1 %i %i\n", v, i*(3) + j+1 ,i*(3) + j +2 );
		}
	}

	fprintf(fpw2,"\nAngles\n\n");
	v=0;
	for(i=0;i<nH2O;i++ ){
		v++;
		fprintf(fpw2,"%i 1 %i %i %i\n", v, i*(3) + 1 ,i*(3) + 2 , i*(3) +3);
	}

	fclose(fpw2);



}