all: moldyn.o posic
-posic:
+posic: moldyn.o $(OBJS)
$(CC) $(CFLAGS) -o $@ $(OBJS) $(LIBS) posic.c
clean:
}
v3_zero(&r);
- v3_add(&r,&r,&o);
count=0;
/* fill up the room */
+ r.x=o.x;
while(r.x<x) {
- r.y=.0;
+ r.y=o.y;
while(r.y<y) {
- r.z=.0;
+ r.z=o.z;
while(r.z<z) {
v3_copy(&(atom[count].r),&r);
+ atom[count].element=1;
count+=1;
for(i=0;i<3;i++) {
v3_add(&n,&r,&basis[i]);
ret=diamond_init(a,b,c,lc,*atom,&origin);
break;
default:
- ret=-1;
printf("unknown lattice type (%02x)\n",type);
+ return -1;
}
/* debug */
printf("ok, there is something wrong ...\n");
printf("calculated -> %d atoms\n",count);
printf("created -> %d atoms\n",ret);
+ return -1;
}
while(count) {
return ret;
}
+int thermal_init(t_atom *atom,int count,double t) {
+
+ /*
+ * - gaussian distribution of velocities
+ * - velocity scaling (E = 3/2 N k T), E: kinetic energy
+ */
+
+ int i;
+ double e,c,v;
+
+ e=.0;
+
+ for(i=0;i<count;i++) {
+ /* x direction */
+ v=gauss_rand();
+ atom[count].v.x=v;
+ e+=0.5*atom[count].mass*v*v;
+ /* y direction */
+ v=gauss_rand();
+ atom[count].v.y=v;
+ e+=0.5*atom[count].mass*v*v;
+ /* z direction */
+ v=gauss_rand();
+ atom[count].v.z=v;
+ e+=0.5*atom[count].mass*v*v;
+ }
+
+ c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN));
+
+ for(i=0;i<count;i++)
+ v3_scale(&(atom[count].v),&(atom[count].v),(1.0/c));
+
+ return 0;
+}
+
/* defines */
+#define K_BOLTZMANN 1.3807E-23;
+
#define FCC 0x01
#define DIAMOND 0x02
/* init */
printf("placing silicon atoms\n");
- count=create_lattice(FCC,Si,M_SI,LC_SI,a,b,c,&si);
+ count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
visual_atoms(&vis,0.0,si,count);
/* script file update */
dprintf(v->fd,"load xyz %s\n",file);
dprintf(v->fd,"spacefill 200\n");
- dprintf(v->fd,"rotate x 11\n");
- dprintf(v->fd,"rotate y 13\n");
+ dprintf(v->fd,"rotate x 100\n");
+ dprintf(v->fd,"rotate y 10\n");
dprintf(v->fd,"set ambient 20\n");
dprintf(v->fd,"set specular on\n");
sprintf(file,"%s-%.15f.ppm",v->fb,time);