y=0.5*dim->y;
z=0.5*dim->z;
- if(a->x>x) a->x-=dim->x;
+ if(a->x>=x) a->x-=dim->x;
else if(-a->x>x) a->x+=dim->x;
- if(a->y>y) a->y-=dim->y;
+ if(a->y>=y) a->y-=dim->y;
else if(-a->y>y) a->y+=dim->y;
- if(a->z>z) a->z-=dim->z;
+ if(a->z>=z) a->z-=dim->z;
else if(-a->z>z) a->z+=dim->z;
return 0;
/* check for visualiziation */
// to be continued ...
- if(!(i%100))
+ if(!(i%1)) {
visual_atoms(moldyn->visual,i*moldyn->tau,
moldyn->atom,moldyn->count);
+ }
}
return 0;
h1*=h2; /* 1/r^14 */
h1*=sig12;
h2*=sig6;
- d=-12.0*h1+6.0*h2;
+ d=12.0*h1-6.0*h2;
d*=eps;
v3_scale(&force,&distance,d);
v3_add(&(atom[j].f),&(atom[j].f),&force);
t=TEMPERATURE;
printf("placing silicon atoms ... ");
- count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
- printf("(%d) ok!\n",count);
+ //count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
+ //printf("(%d) ok!\n",count);
+ count=2;
+ si=malloc(2*sizeof(t_atom));
+ si[0].r.x=2.0;
+ si[0].r.y=0;
+ si[0].r.z=0;
+ si[0].element=Si;
+ si[0].mass=14.0;
+ si[1].r.x=-2.0;
+ si[1].r.y=0;
+ si[1].r.z=0;
+ si[1].element=Si;
+ si[1].mass=14.0;
printf("setting thermal fluctuations\n");
- thermal_init(si,&random,count,t);
-
+ //thermal_init(si,&random,count,t);
+ v3_zero(&(si[0].v));
+ v3_zero(&(si[1].v));
/* check kinetic energy */
md.atom=si;
md.potential=potential_lennard_jones;
md.force=force_lennard_jones;
- md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
+ //md.cutoff_square=((LC_SI/4.0)*(LC_SI/4.0));
+ md.cutoff_square=36.0;
md.pot_params=&lj;
md.integrate=velocity_verlet;
md.time_steps=RUNS;
help=lj.sigma6*lj.sigma6;
lj.sigma6*=help;
lj.sigma12=lj.sigma6*lj.sigma6;
- lj.epsilon=1;
+ lj.epsilon=10000;
u=get_e_pot(&md);
#ifndef POSIC_H
#define POSIC_H
-#define RUNS 15000
-#define TAU 0.000001
+#define RUNS 200
+#define TAU 0.001
#define TEMPERATURE 273.0
-#define LEN_X 1
-#define LEN_Y 1
-#define LEN_Z 1
+#define LEN_X 15
+#define LEN_Y 15
+#define LEN_Z 15
#define R_CUTOFF 20
/* script file update */
dprintf(v->fd,"load xyz %s\n",file);
dprintf(v->fd,"spacefill 200\n");
- dprintf(v->fd,"rotate x 100\n");
- dprintf(v->fd,"rotate y 10\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);