]> hackdaworld.org Git - physik/posic.git/commitdiff
nearly finished init stuff (probs with rand function!)
authorhackbard <hackbard>
Tue, 28 Mar 2006 15:18:41 +0000 (15:18 +0000)
committerhackbard <hackbard>
Tue, 28 Mar 2006 15:18:41 +0000 (15:18 +0000)
Makefile
math/math.c
math/math.h
moldyn.c
moldyn.h
posic.c
posic.h
random/random.c [new file with mode: 0644]
random/random.h [new file with mode: 0644]

index b6e31c03eb7477b32b27a9a476ea703c7240fd3e..ff3fbe3ad5efb511ad2c2d4009cbc211975491c6 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,12 +1,12 @@
 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
index 6d707e93bd7c23b355d8292273349fb37c004d30..9cfaff67296592e49da59f2346e61382ebcf3b3a 100644 (file)
@@ -7,6 +7,8 @@
 
 
 #include <string.h>
+#include <math.h>
+
 #include "math.h"
 
 int v3_add(t_3dvec *sum,t_3dvec *a,t_3dvec *b) {
@@ -64,3 +66,13 @@ int v3_cmp(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)));
+}
+
index df8f6af28c785cc418a9efa1e86e199026bc0233..53e2114641a55f761003750d74af9d37e41eb680 100644 (file)
@@ -25,5 +25,7 @@ int v3_zero(t_3dvec *vec);
 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
index e96cdd4bd66e32632087321d07d12173024bec24..b507d65354dabcf690628bc8bc31846962542fe3 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -9,9 +9,11 @@
 
 #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,
@@ -63,38 +65,83 @@ 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;
+}
+
index 85bbfd73d894115234476408a66f0fd7747c1164..65f0d8da10cfa0753b72ffdeac86dadf7d407917 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -9,6 +9,7 @@
 #define MOLDYN_H
 
 #include "math/math.h"
+#include "random/random.h"
 
 /* datatypes */
 
@@ -23,7 +24,7 @@ typedef struct s_atom {
 
 /* defines */
 
-#define K_BOLTZMANN            1.3807E-23;
+#define K_BOLTZMANN            1.3807E-23
 
 #define FCC                    0x01
 #define DIAMOND                        0x02
@@ -39,5 +40,9 @@ typedef struct s_atom {
 
 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
diff --git a/posic.c b/posic.c
index 0f05b2b6b6ecd2b6bbbe0317efdb9dd7f715ab58..d7a659a34b17466ef1d2d8405677f797236e6f60 100644 (file)
--- a/posic.c
+++ b/posic.c
 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);
 
diff --git a/posic.h b/posic.h
index a26e66830e1cfc2e10b17f3e275ae6c63f87d0da..de0106d84e47c3adfc83ce22c29271ea7d4db842 100644 (file)
--- a/posic.h
+++ b/posic.h
@@ -20,6 +20,8 @@
 #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
diff --git a/random/random.c b/random/random.c
new file mode 100644 (file)
index 0000000..c06b98a
--- /dev/null
@@ -0,0 +1,140 @@
+/*
+ * 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;
+}
+
diff --git a/random/random.h b/random/random.h
new file mode 100644 (file)
index 0000000..8a21276
--- /dev/null
@@ -0,0 +1,48 @@
+/*
+ * 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 */