From 710717c4033bc5b8eb34644914e762a2834ae345 Mon Sep 17 00:00:00 2001 From: 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 + * + */ + +#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.xx) { + while(r.yy) { + while(r.zz) { + v3_copy(&(atom[count].r),&r); + count+=1; + for(i=0;i<3;i++) { + v3_add(&n,&r,&basis[i]); + if((n.xx+o.x)&&(n.yy+o.y)&&(n.zz+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 + * + */ + +#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 + * + */ + + +#include +#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 + * + */ + +#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 + * + */ + +#include "moldyn.h" + +#include +#include + +#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 + * + */ + +#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;iLX) 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;iLX) 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;iLX) 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;iLX) 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 + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include + +#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 + * + */ + +#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