From 710717c4033bc5b8eb34644914e762a2834ae345 Mon Sep 17 00:00:00 2001
From: hackbard <hackbard>
Date: Fri, 24 Mar 2006 19:49:37 +0000
Subject: [PATCH] initial checkin of the new concept approach

---
 Makefile        |  12 ++
 doit            |   2 +-
 init/Makefile   |   7 ++
 init/init.c     |  84 +++++++++++++
 init/init.h     |  20 ++++
 math/Makefile   |   7 ++
 math/math.c     |  66 +++++++++++
 math/math.h     |  29 +++++
 moldyn.c        |  64 ++++++++++
 moldyn.h        |  41 +++++++
 posic.c         | 304 +++++++++++++++++++++---------------------------
 posic.h         |  14 +--
 visual/Makefile |   7 ++
 visual/visual.c |  99 ++++++++++++++++
 visual/visual.h |  25 ++++
 15 files changed, 598 insertions(+), 183 deletions(-)
 create mode 100644 Makefile
 create mode 100644 init/Makefile
 create mode 100644 init/init.c
 create mode 100644 init/init.h
 create mode 100644 math/Makefile
 create mode 100644 math/math.c
 create mode 100644 math/math.h
 create mode 100644 moldyn.c
 create mode 100644 moldyn.h
 create mode 100644 visual/Makefile
 create mode 100644 visual/visual.c
 create mode 100644 visual/visual.h

diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000..a6ba055
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,12 @@
+CC=gcc
+CFLAGS=-Wall
+
+OBJS=init/init.o visual/visual.o math/math.o moldyn.o
+
+all: moldyn.o posic
+
+posic:
+	$(CC) $(CFLAGS) -o $@ $(OBJS) $(LIBS) posic.c
+
+clean:
+	rm -f *.o posic
diff --git a/doit b/doit
index 04ece04..b9bfa91 100755
--- a/doit
+++ b/doit
@@ -4,6 +4,6 @@
 gcc -o posic posic.c -Wall
 ./posic
 ./perms
-rasmol
+rasmol -nodisplay
 ./ppm2avi
 mplayer video/md.avi
diff --git a/init/Makefile b/init/Makefile
new file mode 100644
index 0000000..17a2685
--- /dev/null
+++ b/init/Makefile
@@ -0,0 +1,7 @@
+CC=gcc
+CFLAGS=-Wall
+
+all: init.o ../math/math.o
+
+clean:
+	rm -f *.o 
diff --git a/init/init.c b/init/init.c
new file mode 100644
index 0000000..1f62963
--- /dev/null
+++ b/init/init.c
@@ -0,0 +1,84 @@
+/*
+ * init.c - create initial conditions
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#include "../math/math.h"
+#include "../moldyn.h"
+
+/* fcc lattice init */
+int fcc_init(t_3dvec *dim,double lc,t_atom *atom,t_3dvec *origin) {
+
+	int count;
+	int i,j;
+	t_3dvec o,r,n;
+	t_3dvec basis[3];
+	double help[3];
+
+	if(origin) v3_copy(&o,origin);
+	else v3_zero(&o);
+
+	/* construct the basis */
+	for(i=0;i<3;i++) {
+		for(j=0;j<3;j++) {
+			if(i!=j) help[j]=0.5*lc;
+			else help[j]=.0;
+		}
+		v3_set(&basis[i],help);
+	}
+
+	v3_zero(&r);
+	v3_add(&r,&r,&o);
+	count=0;
+	
+	/* fill up the room */
+	while(r.x<dim->x) {
+		while(r.y<dim->y) {
+			while(r.z<dim->z) {
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+				for(i=0;i<3;i++) {
+					v3_add(&n,&r,&basis[i]);
+					if((n.x<dim->x+o.x)&&(n.y<dim->y+o.y)&&(n.z<dim->z+o.z)) {
+						v3_copy(&(atom[count].r),&n);
+						count+=1;
+					}
+				}
+				r.z+=lc;	
+			}
+			r.y+=lc;
+		}
+		r.x+=lc;
+	}
+
+	/* coordinate transformation */
+	help[0]=dim->x/2.0;
+	help[1]=dim->y/2.0;
+	help[2]=dim->z/2.0;
+	v3_set(&n,help);
+	for(i=0;i<count;i++)
+		v3_sub(&(atom[i].r),&(atom[i].r),&n);
+		
+	return count;
+}
+
+int diamond_init(t_3dvec *dim,double lc,t_atom *atom,t_3dvec *origin) {
+
+	int count;
+	t_3dvec o;
+
+	count=fcc_init(dim,lc,atom,origin);
+
+	o.x=0.25*lc;
+	o.y=0.25*lc;
+	o.z=0.25*lc;
+
+	if(origin) v3_add(&o,&o,origin);
+
+	count+=fcc_init(dim,lc,&atom[count],&o);
+
+
+	return count;
+}
diff --git a/init/init.h b/init/init.h
new file mode 100644
index 0000000..81bfe19
--- /dev/null
+++ b/init/init.h
@@ -0,0 +1,20 @@
+/*
+ * init.h - create initial conditions
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef INIT_H
+#define INIT_H
+
+#include "../moldyn.h"
+
+/* datattypes */
+
+/* function prototypes */
+
+int fcc_init(t_3dvec *dim,double lc,t_atom *atom,t_3dvec *origin);
+int diamond_init(t_3dvec *dim,double lc,t_atom *atom,t_3dvec *origin);
+
+#endif
diff --git a/math/Makefile b/math/Makefile
new file mode 100644
index 0000000..b29746d
--- /dev/null
+++ b/math/Makefile
@@ -0,0 +1,7 @@
+CC=gcc
+CFLAGS=-Wall
+
+all: math.o
+
+clean:
+	rm -f *.o 
diff --git a/math/math.c b/math/math.c
new file mode 100644
index 0000000..6d707e9
--- /dev/null
+++ b/math/math.c
@@ -0,0 +1,66 @@
+/*
+ * math.c - basic vector calculations
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+
+#include <string.h>
+#include "math.h"
+
+int v3_add(t_3dvec *sum,t_3dvec *a,t_3dvec *b) {
+
+	sum->x=a->x+b->x;
+	sum->y=a->y+b->y;
+	sum->z=a->z+b->z;
+
+	return 0;
+}
+
+int v3_sub(t_3dvec *sum,t_3dvec *a,t_3dvec *b) {
+
+	sum->x=a->x-b->x;
+	sum->y=a->y-b->y;
+	sum->z=a->z-b->z;
+
+	return 0;
+}
+
+int v3_scale(t_3dvec *result,t_3dvec *a,double s) {
+
+	result->x=s*a->x;
+	result->y=s*a->y;
+	result->z=s*a->z;
+
+	return 0;
+}
+
+int v3_zero(t_3dvec *vec) {
+
+	vec->x=.0;
+	vec->y=.0;
+	vec->z=.0;
+
+	return 0;
+}
+
+int v3_set(t_3dvec *vec,double *ptr) {
+
+	memcpy(vec,ptr,sizeof(t_3dvec));
+
+	return 0;
+}
+
+int v3_copy(t_3dvec *trg,t_3dvec *src) {
+
+	memcpy(trg,src,sizeof(t_3dvec));
+
+	return 0;
+}
+
+int v3_cmp(t_3dvec *a,t_3dvec *b) {
+
+	return(memcmp(a,b,sizeof(t_3dvec)));
+}
+
diff --git a/math/math.h b/math/math.h
new file mode 100644
index 0000000..df8f6af
--- /dev/null
+++ b/math/math.h
@@ -0,0 +1,29 @@
+/*
+ * math.h - basic vector calculations
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef MATH_H
+#define MATH_H
+
+/* datatypes */
+
+typedef struct s_3dvec {
+	double x;
+	double y;
+	double z;
+} t_3dvec;
+
+/* function prototypes */
+
+int v3_add(t_3dvec *sum,t_3dvec *a,t_3dvec *b);
+int v3_sub(t_3dvec *sum,t_3dvec *a,t_3dvec *b);
+int v3_scale(t_3dvec *result,t_3dvec *a,double s);
+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);
+
+#endif
diff --git a/moldyn.c b/moldyn.c
new file mode 100644
index 0000000..029f7c7
--- /dev/null
+++ b/moldyn.c
@@ -0,0 +1,64 @@
+/*
+ * moldyn.c - molecular dynamics library main file
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#include "moldyn.h"
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "math/math.h"
+#include "init/init.h"
+
+
+int create_lattice(unsigned char type,int element,double mass,double lc,
+                   t_3dvec *dim,t_atom **atom) {
+
+	int count;
+	int ret;
+	t_3dvec origin;
+
+	count=((dim->x/lc)*(dim->y/lc)*(dim->z/lc));
+
+	if(type==FCC) count*=4;
+	if(type==DIAMOND) count*=8;
+
+	*atom=malloc(count*sizeof(t_atom));
+	if(*atom==NULL) {
+		perror("malloc (atoms)");
+		return -1;
+	}
+
+	v3_zero(&origin);
+
+	switch(type) {
+		case FCC:
+			ret=fcc_init(dim,lc,*atom,&origin);
+			break;
+		case DIAMOND:
+			ret=diamond_init(dim,lc,*atom,&origin);
+			break;
+		default:
+			ret=-1;
+			printf("unknown lattice type (%02x)\n",type);
+	}
+
+	/* debug */
+	if(ret!=count) {
+		printf("ok, there is something wrong ...\n");
+		printf("calculated -> %d atoms\n",count);
+		printf("created -> %d atoms\n",ret);
+	}
+
+	while(count) {
+		(*atom)[count-1].element=element;
+		(*atom)[count-1].mass=mass;
+		count-=1;
+	}
+
+	return ret;
+}
+
diff --git a/moldyn.h b/moldyn.h
new file mode 100644
index 0000000..d72ac16
--- /dev/null
+++ b/moldyn.h
@@ -0,0 +1,41 @@
+/*
+ * moldyn.h - molecular dynamics library header file
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef MOLDYN_H
+#define MOLDYN_H
+
+#include "math/math.h"
+
+/* datatypes */
+
+typedef struct s_atom {
+	t_3dvec r;	/* positions */
+	t_3dvec v;	/* velocities */
+	t_3dvec f;	/* forces */
+	int element;	/* number of element in pse */
+	double mass;	/* atom mass */
+} t_atom;
+
+
+/* defines */
+
+#define FCC			0x01
+#define DIAMOND			0x02
+
+#define C			0x06
+#define M_C			6.0
+
+#define Si			0x0e
+#define M_SI			14.0
+#define LC_SI			5.43105
+
+/* function prototypes */
+
+int create_lattice(unsigned char type,int element,double mass,double lc,
+                   t_3dvec *dim,t_atom **atom);
+
+#endif
diff --git a/posic.c b/posic.c
index 8910ddd..6973278 100644
--- a/posic.c
+++ b/posic.c
@@ -5,77 +5,42 @@
  *
  */
  
-#include "posic.h"
+#include "moldyn.h"
+#include "math/math.h"
+#include "init/init.h"
+#include "visual/visual.h"
 
-#define RAND(max) (max*(0.5-(1.0*rand()/RAND_MAX+1)));
+#include "posic.h"
 
 int main(int argc,char **argv) {
 
 	t_atom *si;
-	//t_si *c;
-	int i,j,runs,amount_si;
-	double time;
-	int fd1,fd2;
-	unsigned char xyz[128];
-	unsigned char scr[128];
-	unsigned char ppm[128];
+	t_3dvec dim;
+
+	char fb[32]="saves/fcc_test";
 
-	double tau,tau2,m,m2;
-	double deltax,deltay,deltaz,distance;
-	double deltax2,deltay2,deltaz2,tmp;
-	double lj1,lj2,lj;
+	t_visual vis;
 
-	/* silicon */
-	//amount_si=AMOUNT_SI;
-	amount_si=2;
-	printf("simulating %d silicon atoms\n",amount_si);
-	si=malloc(amount_si*sizeof(t_atom));
-	if(!si) {
-		perror("silicon malloc");
-		return -1;
-	}
-	memset(si,0,amount_si*sizeof(t_atom));
-	m=SI_M; m2=2.0*m;
+	int count;
+
+	dim.x=LEN_X;
+	dim.y=LEN_Y;
+	dim.z=LEN_Z;
+	
+	visual_init(&vis,fb);
 
 	/* init */
 	printf("placing silicon atoms\n");
-	/*
-	for(i=0;i<amount_si;i++) {
-		si[i].x=RAND(LEN_X);
-		si[i].y=RAND(LEN_Y);
-		si[i].z=RAND(LEN_Z);
-		si[i].vx=.0;
-		si[i].vy=.0;
-		si[i].vz=.0;
-		si[i].fx=.0;
-		si[i].fy=.0;
-		si[i].fz=.0;
-	}
-	*/
-	si[0].x=-2.0;
-	si[1].x=+2.0;
-	si[0].y=0;
-	si[0].z=0;
-	si[1].y=0;
-	si[1].z=0;
+	count=create_lattice(FCC,Si,M_SI,LC_SI,&dim,&si);
 
-	/* time */
-	time=.0;
-	tau=TAU;
-	tau2=tau*tau;
+	visual_atoms(&vis,0.0,si,count);
 
-	/* rasmol script file */
-	sprintf(scr,"./saves/si.scr");
-	fd2=open(scr,O_WRONLY|O_CREAT|O_TRUNC);
-	if(fd2<0) {
-		perror("rasmol script file open");
-		return -1;
-	}
+	visual_tini(&vis);
 
-	printf("starting velocity verlet: ");
-	fflush(stdout);
+	//printf("starting velocity verlet: ");
+	//fflush(stdout);
 
-	for(runs=0;runs<RUNS;runs++) {
+	//for(runs=0;runs<RUNS;runs++) {
 
 	/* 
 	 * velocity verlet
@@ -84,118 +49,119 @@ int main(int argc,char **argv) {
 	 * dr/dt|(t+h) = dr/dt|t + h/2m * (F(t) + F(t+h))
 	 *
 	 */
-	for(i=0;i<amount_si;i++) {
-		/* calculation of new positions r(t+h) */
-		si[i].x+=si[i].vx*tau;
-		si[i].y+=si[i].vy*tau;
-		si[i].z+=si[i].vz*tau;
-		si[i].x+=(tau2*si[i].fx/m2);
-		if(si[i].x>LX) si[i].x-=LEN_X;
-		else if(si[i].x<-LX) si[i].x+=LEN_X;
-		si[i].y+=(tau2*si[i].fy/m2);
-		if(si[i].y>LY) si[i].y-=LEN_Y;
-		else if(si[i].y<-LY) si[i].y+=LEN_Y;
-		si[i].z+=(tau2*si[i].fz/m2);
-		if(si[i].z>LZ) si[i].z-=LEN_Z;
-		else if(si[i].z<-LZ) si[i].z+=LEN_Z;
-		/* calculation of velocities v(t+h/2) */
-		si[i].vx+=(tau*si[i].fx/m2);
-		si[i].vy+=(tau*si[i].fy/m2);
-		si[i].vz+=(tau*si[i].fz/m2);
-		/* reset of forces */
-		si[i].fx=.0;
-		si[i].fy=.0;
-		si[i].fz=.0;
-	}
-	for(i=0;i<amount_si;i++) {
-		/* calculation of forces at new positions r(t+h) */
-		for(j=0;j<i;j++) {
-			deltax=si[i].x-si[j].x;
-			if(deltax>LX) deltax-=LEN_X;
-			else if(-deltax>LX) deltax+=LEN_X;
-			deltax2=deltax*deltax;
-			deltay=si[i].y-si[j].y;
-			if(deltay>LY) deltay-=LEN_Y;
-			else if(-deltay>LY) deltay+=LEN_Y;
-			deltay2=deltay*deltay;
-			deltaz=si[i].z-si[j].z;
-			if(deltaz>LZ) deltaz-=LEN_Z;
-			else if(-deltaz>LZ) deltaz+=LEN_Z;
-			deltaz2=deltaz*deltaz;
-			distance=deltax2+deltay2+deltaz2;
-			if(distance<=R2_CUTOFF) {
-				tmp=1.0/distance; // 1/r^2
-				lj1=tmp; // 1/r^2
-				tmp*=tmp; // 1/r^4
-				lj1*=tmp; // 1/r^6
-				tmp*=tmp; // 1/r^8
-				lj2=tmp; // 1/r^8
-				lj1*=tmp; // 1/r^14
-				lj1*=LJ_SIGMA_12;
-				lj2*=LJ_SIGMA_06;
-				lj=-2*lj1+lj2;
-				si[i].fx-=lj*deltax;
-				si[i].fy-=lj*deltay;
-				si[i].fz-=lj*deltaz;
-				si[j].fx+=lj*deltax;
-				si[j].fy+=lj*deltay;
-				si[j].fz+=lj*deltaz;
-			}
-		}
-	}
-	for(i=0;i<amount_si;i++) {
-		/* calculation of new velocities v(t+h) */
-		si[i].vx+=(tau*si[i].fx/m2);
-		si[i].vy+=(tau*si[i].fy/m2);
-		si[i].vz+=(tau*si[i].fz/m2);
-	}
-
-	if(!(runs%100)) {
-
-	/* rasmol script & xyz file */
-	sprintf(xyz,"./saves/si-%.15f.xyz",time);
-	sprintf(ppm,"./video/si-%.15f.ppm",time);
-	fd1=open(xyz,O_WRONLY|O_CREAT|O_TRUNC);
-	if(fd1<0) {
-		perror("rasmol xyz file open");
-		return -1;
-	}
-	dprintf(fd2,"load xyz %s\n",xyz);
-	//dprintf(fd2,"spacefill 200\n");
-	dprintf(fd2,"spacefill\n");
-	dprintf(fd2,"rotate x 11\n");
-	dprintf(fd2,"rotate y 13\n");
-	dprintf(fd2,"set ambient 20\n");
-	dprintf(fd2,"set specular on\n");
-	dprintf(fd2,"write ppm %s\n",ppm);
-	dprintf(fd2,"zap\n");
-	dprintf(fd1,"%d\nsilicon\n",amount_si+9);
-	for(i=0;i<amount_si;i++)
-		dprintf(fd1,"Si %f %f %f %f\n",
-		        si[i].x,si[i].y,si[i].z,time);
-	dprintf(fd1,"H 0.0 0.0 0.0 %f\n",time);
-	dprintf(fd1,"He %f %f %f %f\n",LX,LY,LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",LX,LY,-LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,-LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,-LZ,time);
-	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,-LZ,time);
-	close(fd1);
-
-	}
-
-	/* increase time */
-	time+=tau;
-	printf(".");
-	fflush(stdout);
-
-	}
-
-	printf(" done\n");
-	close(fd2);
-	free(si);
+	//for(i=0;i<amount_si;i++) {
+//		/* calculation of new positions r(t+h) */
+//		si[i].x+=si[i].vx*tau;
+//		si[i].y+=si[i].vy*tau;
+//		si[i].z+=si[i].vz*tau;
+//		si[i].x+=(tau2*si[i].fx/m2);
+//		if(si[i].x>LX) si[i].x-=LEN_X;
+//		else if(si[i].x<-LX) si[i].x+=LEN_X;
+//		si[i].y+=(tau2*si[i].fy/m2);
+//		if(si[i].y>LY) si[i].y-=LEN_Y;
+//		else if(si[i].y<-LY) si[i].y+=LEN_Y;
+//		si[i].z+=(tau2*si[i].fz/m2);
+//		if(si[i].z>LZ) si[i].z-=LEN_Z;
+//		else if(si[i].z<-LZ) si[i].z+=LEN_Z;
+//		/* calculation of velocities v(t+h/2) */
+//		si[i].vx+=(tau*si[i].fx/m2);
+//		si[i].vy+=(tau*si[i].fy/m2);
+//		si[i].vz+=(tau*si[i].fz/m2);
+//		/* reset of forces */
+//		si[i].fx=.0;
+//		si[i].fy=.0;
+//		si[i].fz=.0;
+//	}
+//	for(i=0;i<amount_si;i++) {
+//		/* calculation of forces at new positions r(t+h) */
+//		for(j=0;j<i;j++) {
+//			deltax=si[i].x-si[j].x;
+//			if(deltax>LX) deltax-=LEN_X;
+//			else if(-deltax>LX) deltax+=LEN_X;
+//			deltax2=deltax*deltax;
+//			deltay=si[i].y-si[j].y;
+//			if(deltay>LY) deltay-=LEN_Y;
+//			else if(-deltay>LY) deltay+=LEN_Y;
+//			deltay2=deltay*deltay;
+//			deltaz=si[i].z-si[j].z;
+//			if(deltaz>LZ) deltaz-=LEN_Z;
+//			else if(-deltaz>LZ) deltaz+=LEN_Z;
+//			deltaz2=deltaz*deltaz;
+//			distance=deltax2+deltay2+deltaz2;
+//			if(distance<=R2_CUTOFF) {
+//				tmp=1.0/distance; // 1/r^2
+//				lj1=tmp; // 1/r^2
+//				tmp*=tmp; // 1/r^4
+//				lj1*=tmp; // 1/r^6
+//				tmp*=tmp; // 1/r^8
+//				lj2=tmp; // 1/r^8
+//				lj1*=tmp; // 1/r^14
+//				lj1*=LJ_SIGMA_12;
+//				lj2*=LJ_SIGMA_06;
+//				lj=-2*lj1+lj2;
+//				si[i].fx-=lj*deltax;
+//				si[i].fy-=lj*deltay;
+//				si[i].fz-=lj*deltaz;
+//				si[j].fx+=lj*deltax;
+//				si[j].fy+=lj*deltay;
+//				si[j].fz+=lj*deltaz;
+//			}
+//		}
+//	}
+//	for(i=0;i<amount_si;i++) {
+//		/* calculation of new velocities v(t+h) */
+//		si[i].vx+=(tau*si[i].fx/m2);
+//		si[i].vy+=(tau*si[i].fy/m2);
+//		si[i].vz+=(tau*si[i].fz/m2);
+//	}
+//
+//	if(!(runs%150)) {
+//
+//	/* rasmol script & xyz file */
+//	sprintf(xyz,"./saves/si-%.15f.xyz",time);
+//	sprintf(ppm,"./video/si-%.15f.ppm",time);
+//	fd1=open(xyz,O_WRONLY|O_CREAT|O_TRUNC);
+//	if(fd1<0) {
+//		perror("rasmol xyz file open");
+//		return -1;
+//	}
+//	dprintf(fd2,"load xyz %s\n",xyz);
+//	dprintf(fd2,"spacefill 200\n");
+//	dprintf(fd2,"rotate x 11\n");
+//	dprintf(fd2,"rotate y 13\n");
+//	dprintf(fd2,"set ambient 20\n");
+//	dprintf(fd2,"set specular on\n");
+//	dprintf(fd2,"zoom 400\n");
+//	dprintf(fd2,"write ppm %s\n",ppm);
+//	dprintf(fd2,"zap\n");
+//	dprintf(fd1,"%d\nsilicon\n",amount_si+9);
+//	for(i=0;i<amount_si;i++)
+//		dprintf(fd1,"Si %f %f %f %f\n",
+//		        si[i].x,si[i].y,si[i].z,time);
+//	dprintf(fd1,"H 0.0 0.0 0.0 %f\n",time);
+//	dprintf(fd1,"He %f %f %f %f\n",LX,LY,LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",LX,LY,-LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",-LX,LY,-LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",LX,-LY,-LZ,time);
+//	dprintf(fd1,"He %f %f %f %f\n",-LX,-LY,-LZ,time);
+//	close(fd1);
+//
+//	}
+//
+//	/* increase time */
+//	time+=tau;
+//	printf(".");
+//	fflush(stdout);
+//
+//	}
+//
+//	printf(" done\n");
+//	close(fd2);
+//	free(si);
+//
 
 	return 0;
 }
diff --git a/posic.h b/posic.h
index 755ddc6..1eb8e81 100644
--- a/posic.h
+++ b/posic.h
@@ -24,7 +24,7 @@
 #define LEN_Z 50
 #define LZ (1.0*LEN_Z/2)
 
-#define RUNS 10000
+#define RUNS 15000
 #define TAU 0.001
 
 #define R_CUTOFF 20
@@ -39,16 +39,4 @@
 
 #define AMOUNT_SI ((LEN_X/SI_LC)*(LEN_Y/SI_LC)*(LEN_Z/SI_LC)*2)
 
-typedef struct s_atom {
-	double x;
-	double y;
-	double z;
-	double vx;
-	double vy;
-	double vz;
-	double fx;
-	double fy;
-	double fz;
-} t_atom;
-
 #endif
diff --git a/visual/Makefile b/visual/Makefile
new file mode 100644
index 0000000..b9e1956
--- /dev/null
+++ b/visual/Makefile
@@ -0,0 +1,7 @@
+CC=gcc
+CFLAGS=-Wall
+
+all: visual.o ../math/math.o
+
+clean:
+	rm -f *.o 
diff --git a/visual/visual.c b/visual/visual.c
new file mode 100644
index 0000000..6ad527a
--- /dev/null
+++ b/visual/visual.c
@@ -0,0 +1,99 @@
+/*
+ * visual.c - visualization functions
+ *
+ * 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 <string.h>
+
+#include "visual.h"
+#include "../moldyn.h"
+#include "../math/math.h"
+
+int visual_init(t_visual *v,char *filebase) {
+
+	char file[128+8];
+
+	strncpy(v->fb,filebase,128);
+	memset(file,0,128+8);
+	sprintf(file,"%s.scr",v->fb);
+
+	v->fd=open(file,O_WRONLY);
+	if(v->fd<0) {
+		perror("open visual fd");
+		return -1;
+	}
+
+	return 0;
+}
+
+int visual_tini(t_visual *v) {
+
+	if(v->fd) close(v->fd);
+
+	return 0;
+}
+
+int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
+
+	int i,fd;
+ 	char file[128+64];
+
+	char pse[19][4]={
+		"foo",
+		"H",
+		"He",
+		"Li",
+		"Be",
+		"B",
+		"C",
+		"N",
+		"O",
+		"F",
+		"Ne",
+		"Na",
+		"Mg",
+		"Al",
+		"Si",
+		"P",
+		"S",
+		"Cl",
+		"Ar",
+	};
+
+	sprintf(file,"%s-%.15f.xyz",v->fb,time);
+	fd=open(file,O_WRONLY);
+	if(fd<0) {
+		perror("open visual save file fd");
+		return -1;
+	}
+
+	/* script file update */
+	dprintf(v->fd,"load xyz %s\n",file);
+	dprintf(v->fd,"spacefill 200\n");
+	dprintf(v->fd,"rotate x 11\n");
+	dprintf(v->fd,"rotate y 13\n");
+	dprintf(v->fd,"set ambient 20\n");
+	dprintf(v->fd,"set specular on\n");
+	sprintf(file,"%s-%.15f.ppm",v->fb,time);
+	dprintf(v->fd,"write ppm %s\n",file);
+	dprintf(v->fd,"zap\n");
+
+	/* write the actual data file */
+	dprintf(fd,"Atoms at time %.15f\n",time);
+	dprintf(fd,"%d\n",n);
+	for(i=0;i<n;i++)
+		dprintf(fd,"%s %f %f %f\n",pse[atom[i].element],
+		                           atom[i].r.x,atom[i].r.y,atom[i].r.z);
+	close(fd);
+
+	return 0;
+}
+
diff --git a/visual/visual.h b/visual/visual.h
new file mode 100644
index 0000000..092af84
--- /dev/null
+++ b/visual/visual.h
@@ -0,0 +1,25 @@
+/*
+ * visual.h - visualization header file
+ *                                  
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef VISUAL_H
+#define VISUAL_H
+
+#include "../moldyn.h"
+
+typedef struct s_visual {
+	int fd;			/* rasmol script file descriptor */
+	char fb[128];		/* basename of the save files */
+} t_visual;
+
+
+/* function prototypes */
+
+int visual_init(t_visual *v,char *filebase);
+int visual_tini(t_visual *v);
+int visual_atoms(t_visual *v,double time,t_atom *atom,int n);
+
+#endif
-- 
2.39.2