/* CHANGE LOG April 19,2007 - Initiated this program OBJECTIVE create a LAMMPS data file to put TIP3P water inside a simulation box */ #include #include #include #include #include //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;ipbcl[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;ipbcl[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