$(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
clean:
- rm -f *.o posic
+ rm -f *.o posic */*.o
/* gaussian distribution of velocities */
v3_zero(&p_total);
for(i=0;i<count;i++) {
- sigma=sqrt(2.0*K_BOLTZMANN*t/atom[count].mass);
+ sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
/* x direction */
v=sigma*rand_get_gauss(random);
- atom[count].v.x=v;
- p_total.x+=atom[count].mass*v;
+ atom[i].v.x=v;
+ p_total.x+=atom[i].mass*v;
/* y direction */
v=sigma*rand_get_gauss(random);
- atom[count].v.y=v;
- p_total.x+=atom[count].mass*v;
+ atom[i].v.y=v;
+ p_total.y+=atom[i].mass*v;
/* z direction */
v=sigma*rand_get_gauss(random);
- atom[count].v.z=v;
- p_total.x+=atom[count].mass*v;
+ atom[i].v.z=v;
+ p_total.z+=atom[i].mass*v;
}
/* zero total momentum */
+ v3_scale(&p_total,&p_total,1.0/count);
for(i=0;i<count;i++) {
- v3_scale(&delta,&p_total,1.0/atom[count].mass);
- v3_sub(&(atom[count].v),&(atom[count].v),&delta);
+ v3_scale(&delta,&p_total,1.0/atom[i].mass);
+ v3_sub(&(atom[i].v),&(atom[i].v),&delta);
}
/* velocity scaling */
*/
e=0.0;
for(i=0;i<count;i++)
- e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+ e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN*t));
for(i=0;i<count;i++)
- v3_scale(&(atom[count].v),&(atom[count].v),(1.0/c));
+ v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
return 0;
}
e=0.0;
- for(i=0;i<count;i++)
- e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+ for(i=0;i<count;i++) {
+ e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+ }
return e;
}
+t_3dvec get_total_p(t_atom *atom, int count) {
+
+ t_3dvec p,p_total;
+ int i;
+
+ v3_zero(&p_total);
+ for(i=0;i<count;i++) {
+ v3_scale(&p,&(atom[i].v),atom[i].mass);
+ v3_add(&p_total,&p_total,&p);
+ }
+
+ return p_total;
+}
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
+
+ t_lj_params *params;
+ t_atom *atom;
+ int i,j;
+ int count;
+ t_3dvec distance;
+ double d,help;
+ double u;
+ double eps,sig6,sig12;
+
+ params=ptr;
+ atom=moldyn->atom;
+ count=moldyn->count;
+ eps=params->epsilon;
+ sig6=params->sigma6;
+ sig12=params->sigma12;
+
+ u=0.0;
+ for(i=0;i<count;i++) {
+ for(j=0;j<i;j++) {
+ v3_sub(&distance,&(atom[j].r),&(atom[i].r));
+ d=v3_absolute_square(&distance); /* 1/r^2 */
+ help=d*d; /* 1/r^4 */
+ help*=d; /* 1/r^6 */
+ d=help*help; /* 1/r^12 */
+ u+=eps*(sig12*d-sig6*help);
+ }
+ }
+
+ return u;
+}
+
+
} t_atom;
-/* defines */
+typedef struct s_moldyn {
+ int count;
+ t_atom *atom;
+ double (*potential)(void *ptr);
+ int (*force)(struct s_moldyn *moldyn,void *ptr);
+ unsigned char status;
+} t_moldyn;
-#define K_BOLTZMANN 1.3807E-23
+typedef struct s_lj_params {
+ double sigma6;
+ double sigma12;
+ double epsilon;
+} t_lj_params;
+
+/*
+ * defines
+ */
+
+/* general defines */
+
+#define MOLDYN_STAT_POTENTIAL 0x01
+#define MOLDYN_STAT_FORCE 0x02
+
+/* phsical values */
+
+//#define K_BOLTZMANN 1.3807E-23
+#define K_BOLTZMANN 1.0
#define FCC 0x01
#define DIAMOND 0x02
int thermal_init(t_atom *atom,t_random *random,int count,double t);
int scale_velocity(t_atom *atom,int count,double t);
double get_e_kin(t_atom *atom,int count);
+t_3dvec get_total_p(t_atom *atom,int count);
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr);
#endif
int a,b,c;
double t,e;
+ t_3dvec p;
int count;
char fb[32]="saves/fcc_test";
rand_init(&random,NULL,1);
random.status|=RAND_STAT_VERBOSE;
+ /* testing random numbers */
+ //for(a=0;a<1000000;a++)
+ // printf("%f %f\n",rand_get_gauss(&random),
+ // rand_get_gauss(&random));
+
visual_init(&vis,fb);
a=LEN_X;
t=TEMPERATURE;
- printf("placing silicon atoms\n");
+ printf("placing silicon atoms ... ");
count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
+ printf("(%d) ok!\n",count);
printf("setting thermal fluctuations\n");
thermal_init(si,&random,count,t);
-
/* visualize */
visual_atoms(&vis,0.0,si,count);
printf("kinetic energy: %f\n",e);
printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t);
+ /* check total momentum */
+ p=get_total_p(si,count);
+ printf("total momentum: %f\n",v3_norm(&p));
+
+ /*
+ * let's do the actual md algorithm now
+ *
+ * integration of newtons equations
+ */
+
/* close */
visual_tini(&vis);
#define RUNS 15000
#define TAU 0.001
-#define TEMPERATURE 0.0
+#define TEMPERATURE 273.0
#define SI_M 1
#define SI_LC 5.43105
#define LJ_SIGMA_06 (LJ_SIGMA_02*LJ_SIGMA_02*LJ_SIGMA_02)
#define LJ_SIGMA_12 (LJ_SIGMA_06*LJ_SIGMA_06)
-#define LEN_X 2
-#define LEN_Y 2
-#define LEN_Z 2
+#define LEN_X 1
+#define LEN_Y 1
+#define LEN_Z 1
#define R_CUTOFF 20
#define R2_CUTOFF (R_CUTOFF*R_CUTOFF)
unsigned int rand_get(t_random *random) {
+ int left;
+
if(random->b_ptr==random->buffer+RAND_BUFSIZE) {
if(random->status&RAND_STAT_VERBOSE)
dprintf(random->logfd,
"[random] getting new random numbers\n");
- random->b_ptr=random->buffer;
if(!(random->status&RAND_STAT_UDEV)) {
lseek(random->fd,0,SEEK_SET);
dprintf(random->logfd,
"[random] warning, rereading random file\n");
}
- read(random->fd,random->b_ptr,
- RAND_BUFSIZE*sizeof(unsigned int));
+ left=RAND_BUFSIZE*sizeof(unsigned int);
+ while(left) {
+ left-=read(random->fd,
+ random->buffer+\
+ RAND_BUFSIZE*sizeof(unsigned int)-left,
+ left);
+ }
if(random->status&RAND_STAT_VERBOSE)
dprintf(random->logfd,
"[random] got new random numbers\n");
+ random->b_ptr=random->buffer;
}
return(*(random->b_ptr++));
double rand_get_double(t_random *random) {
- return(1.0*rand_get(random)/(double)(URAND_MAX+1));
+ return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
}
double rand_get_double_linear(t_random *random,double a,double b) {
#define RAND_STAT_UDEV 0x04
#define RAND_STAT_GAUSS 0x08
-#define RAND_BUFSIZE (1024) /* 16 mega byte */
+#define RAND_BUFSIZE (1024) /* 4 k byte */
#define URAND_MAX 0xffffffff
#define URAND_MAX_PLUS_ONE 0x100000000
projects / physik / posic.git / commitdiff