initial checkin of the new concept approach

diff --git a/Makefile b/Makefile
new file mode 100644 (file)
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 (executable)
--- 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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
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 (file)
--- 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 (file)
--- 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 (file)
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 (file)
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 (file)
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