From 3ffe2a08e25fc091b6241885055450009267e2d8 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Tue, 28 Mar 2006 15:18:41 +0000
Subject: [PATCH] nearly finished init stuff (probs with rand function!)

---
 Makefile        |   4 +-
 math/math.c     |  12 +++++
 math/math.h     |   2 +
 moldyn.c        |  69 ++++++++++++++++++++----
 moldyn.h        |   7 ++-
 posic.c         |  36 +++++++++++--
 posic.h         |   2 +
 random/random.c | 140 ++++++++++++++++++++++++++++++++++++++++++++++++
 random/random.h |  48 +++++++++++++++++
 9 files changed, 301 insertions(+), 19 deletions(-)
 create mode 100644 random/random.c
 create mode 100644 random/random.h

diff --git a/Makefile b/Makefile
index b6e31c0..ff3fbe3 100644
--- 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
diff --git a/math/math.c b/math/math.c
index 6d707e9..9cfaff6 100644
--- a/math/math.c
+++ b/math/math.c
@@ -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)));
+}
+
diff --git a/math/math.h b/math/math.h
index df8f6af..53e2114 100644
--- a/math/math.h
+++ b/math/math.h
@@ -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
diff --git a/moldyn.c b/moldyn.c
index e96cdd4..b507d65 100644
--- 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;
+}
+
diff --git a/moldyn.h b/moldyn.h
index 85bbfd7..65f0d8d 100644
--- 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 0f05b2b..d7a659a 100644
--- a/posic.c
+++ b/posic.c
@@ -15,28 +15,54 @@
 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 a26e668..de0106d 100644
--- 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
index 0000000..c06b98a
--- /dev/null
+++ b/random/random.c
@@ -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
index 0000000..8a21276
--- /dev/null
+++ b/random/random.h
@@ -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 */
-- 
2.20.1