2 * posic.c - precipitation process of silicon carbide in silicon
4 * author: Frank Zirkelbach <hackbard@hackdaworld.org>
10 #define RAND(max) (max*(0.5-(1.0*rand()/RAND_MAX+1)));
12 int main(int argc,char **argv) {
16 int i,j,runs,amount_si;
19 unsigned char xyz[128];
20 unsigned char scr[128];
21 unsigned char ppm[128];
24 double deltax,deltay,deltaz,distance;
25 double deltax2,deltay2,deltaz2,tmp;
29 //amount_si=AMOUNT_SI;
31 printf("simulating %d silicon atoms\n",amount_si);
32 si=malloc(amount_si*sizeof(t_atom));
34 perror("silicon malloc");
37 memset(si,0,amount_si*sizeof(t_atom));
41 printf("placing silicon atoms\n");
43 for(i=0;i<amount_si;i++) {
67 /* rasmol script file */
68 sprintf(scr,"./saves/si.scr");
69 fd2=open(scr,O_WRONLY|O_CREAT|O_TRUNC);
71 perror("rasmol script file open");
75 printf("starting velocity verlet: ");
78 for(runs=0;runs<RUNS;runs++) {
83 * r(t+h) = r(t) + h * dr/dt|t + h^2/2m * F(t)
84 * dr/dt|(t+h) = dr/dt|t + h/2m * (F(t) + F(t+h))
87 for(i=0;i<amount_si;i++) {
88 /* calculation of new positions r(t+h) */
89 si[i].x+=si[i].vx*tau;
90 si[i].y+=si[i].vy*tau;
91 si[i].z+=si[i].vz*tau;
92 si[i].x+=(tau2*si[i].fx/m2);
93 if(si[i].x>LX) si[i].x-=LEN_X;
94 else if(si[i].x<-LX) si[i].x+=LEN_X;
95 si[i].y+=(tau2*si[i].fy/m2);
96 if(si[i].y>LY) si[i].y-=LEN_Y;
97 else if(si[i].y<-LY) si[i].y+=LEN_Y;
98 si[i].z+=(tau2*si[i].fz/m2);
99 if(si[i].z>LZ) si[i].z-=LEN_Z;
100 else if(si[i].z<-LZ) si[i].z+=LEN_Z;
101 /* calculation of velocities v(t+h/2) */
102 si[i].vx+=(tau*si[i].fx/m2);
103 si[i].vy+=(tau*si[i].fy/m2);
104 si[i].vz+=(tau*si[i].fz/m2);
105 /* reset of forces */
110 for(i=0;i<amount_si;i++) {
111 /* calculation of forces at new positions r(t+h) */
113 deltax=si[i].x-si[j].x;
114 if(deltax>LX) deltax-=LEN_X;
115 else if(-deltax>LX) deltax+=LEN_X;
116 deltax2=deltax*deltax;
117 deltay=si[i].y-si[j].y;
118 if(deltay>LY) deltay-=LEN_Y;
119 else if(-deltay>LY) deltay+=LEN_Y;
120 deltay2=deltay*deltay;
121 deltaz=si[i].z-si[j].z;
122 if(deltaz>LZ) deltaz-=LEN_Z;
123 else if(-deltaz>LZ) deltaz+=LEN_Z;
124 deltaz2=deltaz*deltaz;
125 distance=deltax2+deltay2+deltaz2;
126 if(distance<=R2_CUTOFF) {
127 tmp=1.0/distance; // 1/r^2
146 for(i=0;i<amount_si;i++) {
147 /* calculation of new velocities v(t+h) */
148 si[i].vx+=(tau*si[i].fx/m2);
149 si[i].vy+=(tau*si[i].fy/m2);
150 si[i].vz+=(tau*si[i].fz/m2);
155 /* rasmol script & xyz file */
156 sprintf(xyz,"./saves/si-%.15f.xyz",time);
157 sprintf(ppm,"./video/si-%.15f.ppm",time);
158 fd1=open(xyz,O_WRONLY|O_CREAT|O_TRUNC);
160 perror("rasmol xyz file open");
163 dprintf(fd2,"load xyz %s\n",xyz);
164 //dprintf(fd2,"spacefill 200\n");
165 dprintf(fd2,"spacefill\n");
166 dprintf(fd2,"rotate x 11\n");
167 dprintf(fd2,"rotate y 13\n");
168 dprintf(fd2,"set ambient 20\n");
169 dprintf(fd2,"set specular on\n");
170 dprintf(fd2,"write ppm %s\n",ppm);
171 dprintf(fd2,"zap\n");
172 dprintf(fd1,"%d\nsilicon\n",amount_si+9);
173 for(i=0;i<amount_si;i++)
174 dprintf(fd1,"Si %f %f %f %f\n",
175 si[i].x,si[i].y,si[i].z,time);
176 dprintf(fd1,"H 0.0 0.0 0.0 %f\n",time);
177 dprintf(fd1,"He %f %f %f %f\n",LX,LY,LZ,time);
178 dprintf(fd1,"He %f %f %f %f\n",-LX,LY,LZ,time);
179 dprintf(fd1,"He %f %f %f %f\n",LX,-LY,LZ,time);
180 dprintf(fd1,"He %f %f %f %f\n",LX,LY,-LZ,time);
181 dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,LZ,time);
182 dprintf(fd1,"He %f %f %f %f\n",-LX,LY,-LZ,time);
183 dprintf(fd1,"He %f %f %f %f\n",LX,-LY,-LZ,time);
184 dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,-LZ,time);