CC=gcc
CFLAGS=-Wall
-OBJS=init/init.o visual/visual.o math/math.o moldyn.o
+OBJS=init/init.o visual/visual.o math/math.o random/random.o moldyn.o
all: moldyn.o posic
posic: moldyn.o $(OBJS)
- $(CC) $(CFLAGS) -o $@ $(OBJS) $(LIBS) posic.c
+ $(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
clean:
rm -f *.o posic
#include <string.h>
+#include <math.h>
+
#include "math.h"
int v3_add(t_3dvec *sum,t_3dvec *a,t_3dvec *b) {
return(memcmp(a,b,sizeof(t_3dvec)));
}
+double v3_absolute_square(t_3dvec *a) {
+
+ return(a->x*a->x+a->y*a->y+a->z*a->z);
+}
+
+double v3_norm(t_3dvec *a) {
+
+ return(sqrt(v3_absolute_square(a)));
+}
+
int v3_set(t_3dvec *vec,double *ptr);
int v3_copy(t_3dvec *trg,t_3dvec *src);
int v3_cmp(t_3dvec *a,t_3dvec *b);
+double v3_absolute_square(t_3dvec *a);
+double v3_norm(t_3dvec *a);
#endif
#include <stdio.h>
#include <stdlib.h>
+#include <math.h>
#include "math/math.h"
#include "init/init.h"
+#include "random/random.h"
int create_lattice(unsigned char type,int element,double mass,double lc,
return ret;
}
-int thermal_init(t_atom *atom,int count,double t) {
+int destroy_lattice(t_atom *atom) {
+
+ if(atom) free(atom);
+
+ return 0;
+}
+
+int thermal_init(t_atom *atom,t_random *random,int count,double t) {
/*
* - gaussian distribution of velocities
+ * - zero total momentum
* - velocity scaling (E = 3/2 N k T), E: kinetic energy
*/
int i;
- double e,c,v;
-
- e=.0;
+ double v,sigma;
+ t_3dvec p_total,delta;
+ /* gaussian distribution of velocities */
+ v3_zero(&p_total);
for(i=0;i<count;i++) {
+ sigma=sqrt(2.0*K_BOLTZMANN*t/atom[count].mass);
/* x direction */
- v=gauss_rand();
+ v=sigma*rand_get_gauss(random);
atom[count].v.x=v;
- e+=0.5*atom[count].mass*v*v;
+ p_total.x+=atom[count].mass*v;
/* y direction */
- v=gauss_rand();
+ v=sigma*rand_get_gauss(random);
atom[count].v.y=v;
- e+=0.5*atom[count].mass*v*v;
+ p_total.x+=atom[count].mass*v;
/* z direction */
- v=gauss_rand();
+ v=sigma*rand_get_gauss(random);
atom[count].v.z=v;
- e+=0.5*atom[count].mass*v*v;
+ p_total.x+=atom[count].mass*v;
+ }
+
+ /* zero total momentum */
+ 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);
}
- c=sqrt((2.0*e)/(3.0*count*K_BOLTZMANN));
+ /* velocity scaling */
+ scale_velocity(atom,count,t);
+ return 0;
+}
+
+int scale_velocity(t_atom *atom,int count,double t) {
+
+ int i;
+ double e,c;
+
+ /*
+ * - velocity scaling (E = 3/2 N k T), E: kinetic energy
+ */
+ e=0.0;
+ for(i=0;i<count;i++)
+ e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].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));
return 0;
}
+double get_e_kin(t_atom *atom,int count) {
+
+ int i;
+ double e;
+
+ e=0.0;
+
+ for(i=0;i<count;i++)
+ e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+
+ return e;
+}
+
#define MOLDYN_H
#include "math/math.h"
+#include "random/random.h"
/* datatypes */
/* defines */
-#define K_BOLTZMANN 1.3807E-23;
+#define K_BOLTZMANN 1.3807E-23
#define FCC 0x01
#define DIAMOND 0x02
int create_lattice(unsigned char type,int element,double mass,double lc,
int a,int b,int c,t_atom **atom);
+int destroy_lattice(t_atom *atom);
+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);
#endif
int main(int argc,char **argv) {
t_atom *si;
+
+ t_visual vis;
+
+ t_random random;
+
int a,b,c;
+ double t,e;
+ int count;
char fb[32]="saves/fcc_test";
- t_visual vis;
+ /* init */
- int count;
+ rand_init(&random,NULL,1);
+ random.status|=RAND_STAT_VERBOSE;
+
+ visual_init(&vis,fb);
a=LEN_X;
b=LEN_Y;
c=LEN_Z;
-
- visual_init(&vis,fb);
- /* init */
+ t=TEMPERATURE;
+
printf("placing silicon atoms\n");
count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
+ printf("setting thermal fluctuations\n");
+ thermal_init(si,&random,count,t);
+
+
+ /* visualize */
+
visual_atoms(&vis,0.0,si,count);
+ /* check kinetic energy */
+
+ e=get_e_kin(si,count);
+ printf("kinetic energy: %f\n",e);
+ printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t);
+
+ /* close */
+
visual_tini(&vis);
+ rand_close(&random);
+
+
//printf("starting velocity verlet: ");
//fflush(stdout);
#define RUNS 15000
#define TAU 0.001
+#define TEMPERATURE 0.0
+
#define SI_M 1
#define SI_LC 5.43105
#define LJ_SIGMA SI_LC
--- /dev/null
+/*
+ * random.c - basic random deviates
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <stdlib.h>
+#include <math.h>
+
+#include "random.h"
+
+int rand_init(t_random *random,char *randomfile,int logfd) {
+
+ random->status=0;
+
+ if(logfd) random->logfd=logfd;
+
+ if(randomfile==NULL) {
+ random->fd=open("/dev/urandom",O_RDONLY);
+ if(random->fd<0) {
+ perror("open urandom");
+ return -1;
+ }
+ random->status|=RAND_STAT_UDEV;
+ } else {
+ random->fd=open(randomfile,O_RDONLY);
+ if(random->fd<0) {
+ perror("open randomfile");
+ return -1;
+ }
+ }
+ random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int));
+ if(random->buffer==NULL) {
+ perror("malloc random buffer");
+ return -1;
+ }
+ random->b_ptr=random->buffer+RAND_BUFSIZE;
+
+ dprintf(random->logfd,"[random] init successfull\n");
+
+ return 0;
+}
+
+
+int rand_close(t_random *random) {
+
+ dprintf(random->logfd,"[random] shutdown\n");
+
+ if(random->buffer) free(random->buffer);
+ if(random->fd) close(random->fd);
+ if(random->logfd) close(random->logfd);
+
+ return 0;
+}
+
+unsigned int rand_get(t_random *random) {
+
+ 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));
+ if(random->status&RAND_STAT_VERBOSE)
+ dprintf(random->logfd,
+ "[random] got new random numbers\n");
+ }
+
+ return(*(random->b_ptr++));
+}
+
+double rand_get_double(t_random *random) {
+
+ return(1.0*rand_get(random)/(double)(URAND_MAX+1));
+}
+
+double rand_get_double_linear(t_random *random,double a,double b) {
+
+ return((-b+sqrt(b*b+2*a*(b+a/2)*rand_get_double(random)))/a);
+}
+
+double rand_get_gauss(t_random *random) {
+
+ double a,b,w;
+
+ if(random->status&RAND_STAT_GAUSS) {
+ random->status&=~RAND_STAT_GAUSS;
+ return random->gauss;
+ }
+
+ w=0;
+ while((w>=1.0)||(w==0.0)) {
+ a=-2.0*rand_get_double(random)+1.0;
+ b=-2.0*rand_get_double(random)+1.0;
+ w=a*a+b*b;
+ }
+
+ w=sqrt(-2.0*log(w)/w);
+
+ random->gauss=b*w;
+ random->status|=RAND_STAT_GAUSS;
+
+ return(a*w);
+}
+
+unsigned int rand_get_int_linear(t_random *random,unsigned int max,
+ double a,double b) {
+
+ return((int)(rand_get_double_linear(random,a,b)*max));
+}
+
+unsigned int rand_get_reject(t_random *random,unsigned int max_x,
+ unsigned int max_y,unsigned int *graph) {
+
+ unsigned int x,y;
+ unsigned char ok;
+
+ ok=0;
+ while(!ok) {
+ x=(max_x*rand_get_double(random));
+ y=(max_y*rand_get_double(random));
+ if(y<=graph[x]) ok=1;
+ }
+
+ return x;
+}
+
--- /dev/null
+/*
+ * random.c - random functions header file
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef RANDOM_H
+#define RANDOM_H
+
+/* datatypes */
+
+typedef struct s_random {
+ int fd;
+ int logfd;
+ unsigned char status;
+ unsigned int *buffer;
+ unsigned int *b_ptr;
+ double gauss;
+} t_random;
+
+/* defines */
+
+#define RAND_STAT_INITIALIZED 0x01
+#define RAND_STAT_VERBOSE 0x02
+#define RAND_STAT_UDEV 0x04
+#define RAND_STAT_GAUSS 0x08
+
+#define RAND_BUFSIZE (1024) /* 16 mega byte */
+
+#define URAND_MAX 0xffffffff
+#define URAND_MAX_PLUS_ONE 0x100000000
+
+/* function prototypes */
+
+int rand_init(t_random *random,char *randomfile,int logfd);
+int rand_close(t_random *random);
+unsigned int rand_get(t_random *random);
+double rand_get_double(t_random *random);
+double rand_get_double_linear(t_random *random,double a,double b);
+double rand_get_gauss(t_random *random);
+unsigned int rand_get_int_linear(t_random *radnom,unsigned int max,
+ double a,double b);
+unsigned int rand_get_reject(t_random *random,unsigned int max_x,
+ unsigned int max_y,unsigned int *graph);
+
+
+#endif /* RANDOM_H */
projects / physik / posic.git / commitdiff