Mc gen.c
From EVOCD
This is the primary polymer generator C code. This file reads in parameters from mc_gen.input and produces a LAMMPS compatible input file for amorphous polyethylene on a FCC lattice.
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
/*
This code generates a series of amorphous polymer chains on
a FCC lattice
*/
#define PI 3.14159265
#define MAX_X 300
#define MAX_Y 300
#define MAX_Z 300
#define MAX_ATOMS 16000000
/* define data structures */
typedef struct atom_t ATOM;
struct atom_t{
int chain; /*chain number*/
int nic; /*number of atom in chain*/
int x[3]; /*lattice occupied by the atom*/
int spot; /*spot or position of atom in lattice*/
int type; /*atom type*/
};
int site[MAX_X][MAX_Y][MAX_Z][4]; /*0 if site is unoccupied 1 if occupied*/
int loc_dens[MAX_X][MAX_Y][MAX_Z][4]; /*local density of site in terms of number of occupied nearest neighbor sites 0-12*/
int snd_dens[MAX_X][MAX_Y][MAX_Z][4]; /*2nd nearest neighbor density 0-6*/
int natom; /*total number of atoms*/
float density;
int c_length;
int n_chains;
int seed;
int n_incs;
float lat;
float mass;
FILE *fp,*fq;
char *inf,*outf,fnames[1][40];
ATOM atom_list[MAX_ATOMS]; /*the array of atoms, the length will be Natom*/
static void read_input(void);
static void gen_pos(void);
static void write_file(void);
main()
{
int i,j,k;
inf = "mc_gen.input";
outf = (char *)malloc(256);
fp = fopen(inf,"r");
if(fp==NULL) printf("input file not found\n");
printf("enter read _input\n");
read_input();
printf("exit read _input\n");
printf("exit read _input\n");
srand(seed);
printf("enter gen_pos\n");
gen_pos();
printf("enter write_file\n");
fq = fopen(outf,"w");
if(fq==NULL) printf("output file not found\n");
write_file();
printf("exit write_file\n");
fclose(fq);
}
static void read_input(void)
{
int i,j,k,m,u,v;
int ii,jj,kk,mm;
char *keyw,buf[256];
/* Read output file name */
fgets(fnames[0],sizeof(fnames[0]),fp);
/* Read approximate density */
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%f",&density);
/* Read chain length*/
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%d",&c_length);
/* Read number of chains*/
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%d",&n_chains);
/* Read lattice size*/
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%f",&lat);
lat=lat/sqrt(2.0)*2.0;
/* Read random number seed*/
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%d",&seed);
/* Read random number seed*/
fgets(buf,sizeof(buf),fp);
fgets(buf,sizeof(buf),fp);
sscanf(buf,"%f",&mass);
printf("den %f c_leng %d n_ch %d l_s %f seed %d \n",density,c_length,n_chains,lat,seed);
fclose(fp);
printf("create names \n");
printf("fname is %s \n",fnames[0]);
strncpy(outf,fnames[0],strlen(fnames[0])-1);
printf("output file is is %s \n",outf);
}
static void gen_pos(void)
{
int i,j,k,m,u,v;
int ii,jj,kk,mm,nn,n;
int n_oc_sites,t_sites;
int dx[3],dy[3],dz[3];
int st_ld[12][4];
int snd_ld[6][4];
int snd_count,nn_count,total_dens,point;
int prob_low,prob_high;
int reset;
double sqrt_2_2;
double angle,va_mag,vb_mag,d_prod;
double a[3],aa[3],va[3],vb[3],vc[3];
double neigh[12][4],v_store[12][3],mag_store[12];
char buf[256];
/* Determine the edge lengths in lattice units from the density */
sqrt_2_2=sqrt(2)/2.0;
n_oc_sites=c_length*n_chains;
t_sites=(int)(n_oc_sites/density);
n_incs=(int)(cbrt(t_sites/4.0))+1.0;
printf("%f %d %d %d \n",sqrt_2_2,n_oc_sites,t_sites,n_incs);
/* Initialize the sites as empty */
natom=0;
for(i=0;i<n_incs;i++)
{
for(j=0;j<n_incs;j++)
{
for(k=0;k<n_incs;k++)
{
for(m=0;m<4;m++)
{
site[i][j][k][m]=0;
loc_dens[i][j][k][m]=0;
snd_dens[i][j][k][m]=0;
}
}
}
}
for(u=0;u<n_chains;u++)
{
atom_list[natom].chain=u+1;
atom_list[natom].nic=0;
i=(int)n_incs*((double)rand()/(double)RAND_MAX);
j=(int)n_incs*((double)rand()/(double)RAND_MAX);
k=(int)n_incs*((double)rand()/(double)RAND_MAX);
m=(int)4*((double)rand()/(double)RAND_MAX);
if(site[i][j][k][m]==1)
{
u--;
}
else
{
printf("n_incs %d i %d j %d k %d m %d\n",n_incs,i,j,k,m);
atom_list[natom].x[0]=i;
atom_list[natom].x[1]=j;
atom_list[natom].x[2]=k;
atom_list[natom].spot=m;
atom_list[natom].type=1;
site[i][j][k][m]=1;
natom++;
for(v=1;v<c_length;v++)
{
printf("looking for neighbor %d \n",v);
nn=0;
nn_count=0;
snd_count=0;
angle=0;
dx[1]=i;
dy[1]=j;
dz[1]=k;
dx[2]=i+1;
dy[2]=j+1;
dz[2]=k+1;
dx[0]=i-1;
dy[0]=j-1;
dz[0]=k-1;
if(dx[2]==n_incs)
{
dx[2]=0;
}
if(dy[2]==n_incs)
{
dy[2]=0;
}
if(dz[2]==n_incs)
{
dz[2]=0;
}
if(dx[0]<0)
{
dx[0]=n_incs-1;
}
if(dy[0]<0)
{
dy[0]=n_incs-1;
}
if(dz[0]<0)
{
dz[0]=n_incs-1;
}
if(m==0)
{
a[0]=0;
a[1]=0;
a[2]=0;
}
else if(m==1)
{
a[0]=0.5;
a[1]=0.5;
a[2]=0;
}
else if(m==2)
{
a[0]=0.5;
a[1]=0;
a[2]=0.5;
}
else
{
a[0]=0;
a[1]=0.5;
a[2]=0.5;
}
for(ii=0;ii<3;ii++)
{
for(jj=0;jj<3;jj++)
{
for(kk=0;kk<3;kk++)
{
for(mm=0;mm<4;mm++)
{
if(mm==0)
{
aa[0]=0;
aa[1]=0;
aa[2]=0;
}
else if(mm==1)
{
aa[0]=0.5;
aa[1]=0.5;
aa[2]=0;
}
else if(mm==2)
{
aa[0]=0.5;
aa[1]=0;
aa[2]=0.5;
}
else if(mm==3)
{
aa[0]=0;
aa[1]=0.5;
aa[2]=0.5;
}
aa[0]=aa[0]+ii-1;
aa[1]=aa[1]+jj-1;
aa[2]=aa[2]+kk-1;
va[0]=(double)aa[0]-(double)a[0];
va[1]=(double)aa[1]-(double)a[1];
va[2]=(double)aa[2]-(double)a[2];
va_mag=sqrt(va[0]*va[0]+va[1]*va[1]+va[2]*va[2]);
if(site[dx[ii]][dy[jj]][dz[kk]][mm]==0||(v==c_length-1&&site[dx[ii]][dy[jj]][dz[kk]][mm]==2))
{
if(va_mag<=sqrt_2_2+.01)
{
loc_dens[dx[ii]][dy[jj]][dz[kk]][mm]++;
st_ld[nn_count][0]=dx[ii];
st_ld[nn_count][1]=dy[jj];
st_ld[nn_count][2]=dz[kk];
st_ld[nn_count][3]=mm;
nn_count++;
}
/*
else if(va_mag<=1.01)
{
snd_dens[dx[ii]][dy[jj]][dz[kk]][mm]++;
snd_ld[snd_count][0]=dx[ii];
snd_ld[snd_count][1]=dy[jj];
snd_ld[snd_count][2]=dz[kk];
snd_ld[snd_count][3]=mm;
snd_count++;
}
*/
if(va_mag<sqrt_2_2)
{
printf("creating atoms too close to each other\n");
exit(1);
}
if(va_mag<=sqrt_2_2+.01)
{
if(v!=1)
{
d_prod=va[0]*vb[0]+va[1]*vb[1]+va[2]*vb[2];
angle=acos(d_prod/(va_mag*vb_mag))*180.0/PI;
}
if((angle<-60.0&&angle>-150.0)||(angle>60.0&&angle<150.0)||v==1)
{
neigh[nn][0]=dx[ii];
neigh[nn][1]=dy[jj];
neigh[nn][2]=dz[kk];
neigh[nn][3]=mm;
v_store[nn][0]=va[0];
v_store[nn][1]=va[1];
v_store[nn][2]=va[2];
mag_store[nn]=va_mag;
nn++;
}
}
}
}
}
}
}
if(nn==0)
{
if(reset==1)
{
printf("has already been reset and still can not find a neighbor\n");
printf("chain %d atom %d \n",u,v);
exit(1);
}
printf("there is not a free neighbor site and must reset\n");
reset=1;
site[i][j][k][m]=2;
natom--;
i=atom_list[natom].x[0];
j=atom_list[natom].x[1];
k=atom_list[natom].x[2];
m=atom_list[natom].spot;
for(ii=0;ii<12;ii++)
{
loc_dens[st_ld[ii][0]][st_ld[ii][1]][st_ld[ii][2]][st_ld[ii][3]]--;
}
/*
for(ii=0;ii<6;ii++)
{
snd_dens[snd_ld[ii][0]][snd_ld[ii][1]][snd_ld[ii][2]][snd_ld[ii][3]]--;
}
*/
vb[0]=vc[0];
vb[1]=vc[1];
vb[2]=vc[2];
v=v-2;
}
else
{
total_dens=0;
for(ii=0;ii<nn;ii++)
{
total_dens=total_dens+(12-loc_dens[st_ld[ii][0]][st_ld[ii][1]][st_ld[ii][2]][st_ld[ii][3]]);
/*
total_dens=total_dens+(6-snd_dens[st_ld[ii][0]][st_ld[ii][1]][st_ld[ii][2]][st_ld[ii][3]]);
*/
}
point=total_dens*((double)rand()/(double)RAND_MAX);
prob_low=0;
prob_high=(12-loc_dens[st_ld[0][0]][st_ld[0][1]][st_ld[0][2]][st_ld[0][3]]);
/*
prob_high=prob_high+(6-snd_dens[st_ld[0][0]][st_ld[0][1]][st_ld[0][2]][st_ld[0][3]]);
*/
printf("span %d prob_high %d prob_low %d point %d\n",prob_high-prob_low,prob_high,prob_low,point);
printf("0 %d 1 %d 2 %d 3 %d \n",st_ld[0][0],st_ld[0][1],st_ld[0][2],st_ld[0][3]);
printf("site %d\n",site[st_ld[0][0]][st_ld[0][1]][st_ld[0][2]][st_ld[0][3]]);
for(ii=0;ii<nn;ii++)
{
if(point>=prob_low&&point<prob_high)
{
i=neigh[ii][0];
j=neigh[ii][1];
k=neigh[ii][2];
m=neigh[ii][3];
atom_list[natom].x[0]=i;
atom_list[natom].x[1]=j;
atom_list[natom].x[2]=k;
atom_list[natom].spot=m;
atom_list[natom].type=1;
atom_list[natom].chain=u+1;
site[i][j][k][m]=1;
natom++;
/*store the previous vector*/
vc[0]=vb[0];
vc[1]=vb[1];
vc[2]=vb[2];
vb[0]=-v_store[ii][0];
vb[1]=-v_store[ii][1];
vb[2]=-v_store[ii][2];
vb_mag=mag_store[ii];
ii=nn;
reset=0;
}
else
{
if(ii==nn-1)
{
if(reset==1)
{
printf("there is not a free neighbor site and must reset\n");
printf("has already been reset and still can not find a neighbor\n");
printf("chain %d atom %d \n",u,v);
exit(1);
}
printf("resetting\n");
reset=1;
site[i][j][k][m]=2;
natom--;
i=atom_list[natom].x[0];
j=atom_list[natom].x[1];
k=atom_list[natom].x[2];
m=atom_list[natom].spot;
for(ii=0;ii<12;ii++)
{
loc_dens[st_ld[ii][0]][st_ld[ii][1]][st_ld[ii][2]][st_ld[ii][3]]--;
}
/*
for(ii=0;ii<6;ii++)
{
snd_dens[snd_ld[ii][0]][snd_ld[ii][1]][snd_ld[ii][2]][snd_ld[ii][3]]--;
}
*/
vb[0]=vc[0];
vb[1]=vc[1];
vb[2]=vc[2];
ii=nn;
v=v-2;
}
else
{
prob_low=prob_high;
prob_high=prob_high+(12-loc_dens[st_ld[ii+1][0]][st_ld[ii+1][1]][st_ld[ii+1][2]][st_ld[ii+1][3]]);
/*
prob_high=prob_high+(6-snd_dens[st_ld[ii+1][0]][st_ld[ii+1][1]][st_ld[ii+1][2]][st_ld[ii+1][3]]);
*/
}
}
}
}
}
}
}
}
static void write_file(void)
{
int i,j,k;
int atom_1,atom_2,atom_3,atom_4;
int count;
double a[3],xx[3];
/*
fprintf(fq,"TITLE=\"DATA FILE\"\n");
fprintf(fq,"VARIABLES=\"ID\" \"CHAIN\" \"X\" \"Y\" \"Z\"\n");
fprintf(fq,"ZONE F=POINT\n");
*/
printf("entered\n");
fprintf(fq,"# Model for PE\n");
fprintf(fq,"\n");
fprintf(fq,"%10d atoms\n",natom);
fprintf(fq,"%10d bonds\n",natom-n_chains);
fprintf(fq,"%10d angles\n",natom-2*n_chains);
fprintf(fq,"%10d dihedrals\n",natom-3*n_chains);
fprintf(fq,"\n");
fprintf(fq,"%10d atom types\n",1);
fprintf(fq,"%10d bond types\n",1);
fprintf(fq,"%10d angle types\n",1);
fprintf(fq,"%10d dihedral types\n",1);
fprintf(fq,"\n");
fprintf(fq,"%10.4f%10.4f xlo xhi\n",0.0,(double)n_incs*lat);
fprintf(fq,"%10.4f%10.4f ylo yhi\n",0.0,(double)n_incs*lat);
fprintf(fq,"%10.4f%10.4f zlo zhi\n",0.0,(double)n_incs*lat);
fprintf(fq,"\n");
fprintf(fq,"Masses\n");
fprintf(fq,"\n");
fprintf(fq,"%10d %14.2f\n",1,mass);
fprintf(fq,"\n");
fprintf(fq,"Atoms\n");
fprintf(fq,"\n");
for(i=0;i<natom;i++)
{
xx[0]=(double)atom_list[i].x[0];
xx[1]=(double)atom_list[i].x[1];
xx[2]=(double)atom_list[i].x[2];
if(atom_list[i].spot==0)
{
a[0]=0;
a[1]=0;
a[2]=0;
}
else if(atom_list[i].spot==1)
{
a[0]=.5;
a[1]=.5;
a[2]=0;
}
else if(atom_list[i].spot==2)
{
a[0]=.5;
a[1]=0;
a[2]=.5;
}
else if(atom_list[i].spot==3)
{
a[0]=0;
a[1]=.5;
a[2]=.5;
}
for(j=0;j<3;j++)
{
xx[j]=lat*(xx[j]+a[j]);
}
fprintf(fq,"%10d%10d%10d%10.4f%10.4f%10.4f\n",i+1,atom_list[i].chain,atom_list[i].type,xx[0],xx[1],xx[2]);
}
fprintf(fq,"\n");
fprintf(fq,"Bonds \n");
fprintf(fq,"\n");
count=1;
atom_1=1;
atom_2=2;
for(i=0;i<n_chains;i++)
{
for(j=0;j<c_length-1;j++)
{
fprintf(fq,"%10d%10d%10d%10d\n",count,1,atom_1,atom_2);
count++;
atom_1++;
atom_2++;
}
atom_1++;
atom_2++;
}
fprintf(fq,"\n");
fprintf(fq,"Angles \n");
fprintf(fq,"\n");
count=1;
atom_1=1;
atom_2=2;
atom_3=3;
for(i=0;i<n_chains;i++)
{
for(j=0;j<c_length-2;j++)
{
fprintf(fq,"%10d%10d%10d%10d%10d\n",count,1,atom_1,atom_2,atom_3);
count++;
atom_1++;
atom_2++;
atom_3++;
}
atom_1=atom_1+2;
atom_2=atom_2+2;
atom_3=atom_3+2;
}
fprintf(fq,"\n");
fprintf(fq,"Dihedrals \n");
fprintf(fq,"\n");
count=1;
atom_1=1;
atom_2=2;
atom_3=3;
atom_4=4;
for(i=0;i<n_chains;i++)
{
for(j=0;j<c_length-3;j++)
{
fprintf(fq,"%10d%10d%10d%10d%10d%10d\n",count,1,atom_1,atom_2,atom_3,atom_4);
count++;
atom_1++;
atom_2++;
atom_3++;
atom_4++;
}
atom_1=atom_1+3;
atom_2=atom_2+3;
atom_3=atom_3+3;
atom_4=atom_4+3;
}
}
|
