X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=posic.c;h=6b92a5e51b0428de417e34392189a6aa246abffb;hp=006fd99a298211458df528cb3986d3765a138f7a;hb=HEAD;hpb=9ae5e8af9192b47fba3743ab7f4ef6d0c34d507d diff --git a/posic.c b/posic.c index 006fd99..6b92a5e 100644 --- a/posic.c +++ b/posic.c @@ -1,153 +1,41 @@ /* * posic.c - precipitation process of silicon carbide in silicon * - * author: Frank Zirkelbach + * author: Frank Zirkelbach * */ - + +/* main include file */ + #include "posic.h" -#define RAND(max) (max*(0.5-(1.0*rand()/RAND_MAX+1))); +/* functions */ + + + +/* main code */ + +int parse_config_file() { + + return 0; +} int main(int argc,char **argv) { - t_atom *si; - //t_si *c; - int i,j,runs,amount_si; - double time; - int fd; - - double tau,tau2,m,m2; - double deltax,deltay,deltaz,distance; - double deltax2,deltay2,deltaz2,tmp; - double lj1,lj2,lj; - - /* silicon */ - amount_si=AMOUNT_SI; - printf("simulating %d silicon atoms\n",amount_si); - si=malloc(amount_si*sizeof(t_atom)); - if(!si) { - perror("silicon malloc"); - return -1; - } - memset(si,0,amount_si*sizeof(t_atom)); - m=SI_M; m2=2.0*m; - - /* init */ - printf("placing silicon atoms\n"); - for(i=0;iLX) si[i].x-=LEN_X; - else if(si[i].x<-LX) si[i].x+=LEN_X; - si[i].y+=(tau2*si[i].fy/m2); - if(si[i].y>LY) si[i].y-=LEN_Y; - else if(si[i].y<-LY) si[i].y+=LEN_Y; - si[i].z+=(tau2*si[i].fz/m2); - if(si[i].z>LZ) si[i].z-=LEN_Z; - else if(si[i].z<-LZ) si[i].z+=LEN_Z; - /* calculation of velocities v(t+h/2) */ - si[i].vx+=(tau*si[i].fx/m2); - si[i].vy+=(tau*si[i].fy/m2); - si[i].vz+=(tau*si[i].fz/m2); - } - for(i=0;iLX) deltax-=LEN_X; - else if(-deltax>LX) deltax+=LEN_X; - deltax2=deltax*deltax; - deltay=si[i].y-si[j].y; - if(deltay>LY) deltay-=LEN_Y; - else if(-deltay>LY) deltay+=LEN_Y; - deltay2=deltay*deltay; - deltaz=si[i].z-si[j].z; - if(deltaz>LZ) deltaz-=LEN_Z; - else if(-deltaz>LZ) deltaz+=LEN_Z; - deltaz2=deltaz*deltaz; - distance=deltax2+deltay2+deltaz2; - if(distance<=R2_CUTOFF) { - tmp=1.0/distance; // 1/r^2 - lj1=tmp; // 1/r^2 - tmp*=tmp; // 1/r^4 - lj1*=tmp; // 1/r^6 - tmp*=tmp; // 1/r^8 - lj2=tmp; // 1/r^8 - lj1*=tmp; // 1/r^14 - lj1*=LJ_SIGMA_12; - lj2*=LJ_SIGMA_06; - lj=-2*lj1+lj2; - si[i].fx=lj*deltax; - si[i].fy=lj*deltay; - si[i].fz=lj*deltaz; - si[i].fx=-lj*deltax; - si[i].fy=-lj*deltay; - si[i].fz=-lj*deltaz; - } - } - /* calculation of new velocities v(t+h) */ - si[i].vx+=(tau*si[i].fx/m2); - si[i].vy+=(tau*si[i].fy/m2); - si[i].vz+=(tau*si[i].fz/m2); - } - - time+=tau; - - /* print out positions in rasmol format */ - dprintf(fd,"%d\nTime %f\n",amount_si,time); - for(i=0;i