started potentials

diff --git a/Makefile b/Makefile
index ff3fbe3ad5efb511ad2c2d4009cbc211975491c6..854b5e6f527790975ae11fb917e7faf3a75eecd2 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -9,4 +9,4 @@ posic: moldyn.o $(OBJS)
        $(CC) $(CFLAGS) -lm -o $@ $(OBJS) $(LIBS) posic.c
 
 clean:
-       rm -f *.o posic
+       rm -f *.o posic */*.o

diff --git a/moldyn.c b/moldyn.c
index b507d65354dabcf690628bc8bc31846962542fe3..ecba8baffc40a07c641eee9bfed8838c7c8116cc 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -87,25 +87,26 @@ int thermal_init(t_atom *atom,t_random *random,int count,double t) {
        /* gaussian distribution of velocities */
        v3_zero(&p_total);
        for(i=0;i<count;i++) {
-               sigma=sqrt(2.0*K_BOLTZMANN*t/atom[count].mass);
+               sigma=sqrt(2.0*K_BOLTZMANN*t/atom[i].mass);
                /* x direction */
                v=sigma*rand_get_gauss(random);
-               atom[count].v.x=v;
-               p_total.x+=atom[count].mass*v;
+               atom[i].v.x=v;
+               p_total.x+=atom[i].mass*v;
                /* y direction */
                v=sigma*rand_get_gauss(random);
-               atom[count].v.y=v;
-               p_total.x+=atom[count].mass*v;
+               atom[i].v.y=v;
+               p_total.y+=atom[i].mass*v;
                /* z direction */
                v=sigma*rand_get_gauss(random);
-               atom[count].v.z=v;
-               p_total.x+=atom[count].mass*v;
+               atom[i].v.z=v;
+               p_total.z+=atom[i].mass*v;
        }
 
        /* zero total momentum */
+       v3_scale(&p_total,&p_total,1.0/count);
        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);
+               v3_scale(&delta,&p_total,1.0/atom[i].mass);
+               v3_sub(&(atom[i].v),&(atom[i].v),&delta);
        }
 
        /* velocity scaling */
@@ -124,10 +125,10 @@ int scale_velocity(t_atom *atom,int count,double t) {
         */
        e=0.0;
        for(i=0;i<count;i++)
-               e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+               e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].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));
+               v3_scale(&(atom[i].v),&(atom[i].v),(1.0/c));
 
        return 0;
 }
@@ -139,9 +140,58 @@ double get_e_kin(t_atom *atom,int count) {
 
        e=0.0;
 
-       for(i=0;i<count;i++)
-               e+=0.5*atom[count].mass*v3_absolute_square(&(atom[count].v));
+       for(i=0;i<count;i++) {
+               e+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+       }
 
        return e;
 }
 
+t_3dvec get_total_p(t_atom *atom, int count) {
+
+       t_3dvec p,p_total;
+       int i;
+
+       v3_zero(&p_total);
+       for(i=0;i<count;i++) {
+               v3_scale(&p,&(atom[i].v),atom[i].mass);
+               v3_add(&p_total,&p_total,&p);
+       }
+
+       return p_total;
+}
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
+
+       t_lj_params *params;
+       t_atom *atom;
+       int i,j;
+       int count;
+       t_3dvec distance;
+       double d,help;
+       double u;
+       double eps,sig6,sig12;
+
+       params=ptr;
+       atom=moldyn->atom;
+       count=moldyn->count;
+       eps=params->epsilon;
+       sig6=params->sigma6;
+       sig12=params->sigma12;
+
+       u=0.0;
+       for(i=0;i<count;i++) {
+               for(j=0;j<i;j++) {
+                       v3_sub(&distance,&(atom[j].r),&(atom[i].r));
+                       d=v3_absolute_square(&distance);        /* 1/r^2 */
+                       help=d*d;                               /* 1/r^4 */
+                       help*=d;                                /* 1/r^6 */
+                       d=help*help;                            /* 1/r^12 */
+                       u+=eps*(sig12*d-sig6*help);
+               }
+       }
+       
+       return u;
+}
+
+

diff --git a/moldyn.h b/moldyn.h
index 65f0d8da10cfa0753b72ffdeac86dadf7d407917..30c4bf89527b6d023534ca02f2df1845711f70ba 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -22,9 +22,33 @@ typedef struct s_atom {
 } t_atom;
 
 
-/* defines */
+typedef struct s_moldyn {
+       int count;
+       t_atom *atom;
+       double (*potential)(void *ptr);
+       int (*force)(struct s_moldyn *moldyn,void *ptr);
+       unsigned char status;
+} t_moldyn;
 
-#define K_BOLTZMANN            1.3807E-23
+typedef struct s_lj_params {
+       double sigma6;
+       double sigma12;
+       double epsilon;
+} t_lj_params;
+
+/*
+ *  defines
+ */
+
+/* general defines */
+
+#define MOLDYN_STAT_POTENTIAL          0x01
+#define MOLDYN_STAT_FORCE              0x02
+
+/* phsical values */
+
+//#define K_BOLTZMANN          1.3807E-23
+#define K_BOLTZMANN            1.0
 
 #define FCC                    0x01
 #define DIAMOND                        0x02
@@ -44,5 +68,8 @@ 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);
+t_3dvec get_total_p(t_atom *atom,int count);
+
+double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr);
 
 #endif

diff --git a/posic.c b/posic.c
index d7a659a34b17466ef1d2d8405677f797236e6f60..4eb2d043e4be7972848ef059b4f69ad39204ce6d 100644 (file)
--- a/posic.c
+++ b/posic.c
@@ -22,6 +22,7 @@ int main(int argc,char **argv) {
 
        int a,b,c;
        double t,e;
+       t_3dvec p;
        int count;
 
        char fb[32]="saves/fcc_test";
@@ -31,6 +32,11 @@ int main(int argc,char **argv) {
        rand_init(&random,NULL,1);
        random.status|=RAND_STAT_VERBOSE;
 
+       /* testing random numbers */
+       //for(a=0;a<1000000;a++)
+       //      printf("%f %f\n",rand_get_gauss(&random),
+       //                       rand_get_gauss(&random));
+
        visual_init(&vis,fb);
 
        a=LEN_X;
@@ -39,13 +45,13 @@ int main(int argc,char **argv) {
 
        t=TEMPERATURE;
 
-       printf("placing silicon atoms\n");
+       printf("placing silicon atoms ... ");
        count=create_lattice(DIAMOND,Si,M_SI,LC_SI,a,b,c,&si);
+       printf("(%d) ok!\n",count);
 
        printf("setting thermal fluctuations\n");
        thermal_init(si,&random,count,t);
 
-
        /* visualize */
 
        visual_atoms(&vis,0.0,si,count);
@@ -56,6 +62,16 @@ int main(int argc,char **argv) {
        printf("kinetic energy: %f\n",e);
        printf("3/2 N k T = %f\n",1.5*count*K_BOLTZMANN*t);
 
+       /* check total momentum */
+       p=get_total_p(si,count);
+       printf("total momentum: %f\n",v3_norm(&p));
+
+       /*
+        * let's do the actual md algorithm now
+        *
+        * integration of newtons equations
+        */
+
        /* close */
 
        visual_tini(&vis);

diff --git a/posic.h b/posic.h
index de0106d84e47c3adfc83ce22c29271ea7d4db842..8f306ec913bd7c9225b8eff72384e6669c44f5ff 100644 (file)
--- a/posic.h
+++ b/posic.h
@@ -20,7 +20,7 @@
 #define RUNS 15000
 #define TAU 0.001
 
-#define TEMPERATURE 0.0
+#define TEMPERATURE 273.0 
 
 #define SI_M 1
 #define SI_LC 5.43105
@@ -29,9 +29,9 @@
 #define LJ_SIGMA_06 (LJ_SIGMA_02*LJ_SIGMA_02*LJ_SIGMA_02)
 #define LJ_SIGMA_12 (LJ_SIGMA_06*LJ_SIGMA_06)
 
-#define LEN_X 2
-#define LEN_Y 2
-#define LEN_Z 2
+#define LEN_X 1
+#define LEN_Y 1
+#define LEN_Z 1
 
 #define R_CUTOFF 20
 #define R2_CUTOFF (R_CUTOFF*R_CUTOFF)

diff --git a/random/random.c b/random/random.c
index c06b98ad11de103e73f466d2f7b7e6241f1a3e38..ddbd7b5effbd86f02dd15f9f78d54d1a0a81ebe0 100644 (file)
--- a/random/random.c
+++ b/random/random.c
@@ -62,21 +62,28 @@ int rand_close(t_random *random) {
 
 unsigned int rand_get(t_random *random) {
 
+       int left;
+
        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));
+               left=RAND_BUFSIZE*sizeof(unsigned int);
+               while(left) {
+                       left-=read(random->fd,
+                                  random->buffer+\
+                                  RAND_BUFSIZE*sizeof(unsigned int)-left,
+                                  left);
+               }
                if(random->status&RAND_STAT_VERBOSE)
                        dprintf(random->logfd,
                                "[random] got new random numbers\n");
+               random->b_ptr=random->buffer;
        }
 
        return(*(random->b_ptr++));
@@ -84,7 +91,7 @@ unsigned int rand_get(t_random *random) {
 
 double rand_get_double(t_random *random) {
 
-       return(1.0*rand_get(random)/(double)(URAND_MAX+1));
+       return(1.0*rand_get(random)/((long long unsigned int)URAND_MAX+1));
 }
 
 double rand_get_double_linear(t_random *random,double a,double b) {

diff --git a/random/random.h b/random/random.h
index 8a212767885730f10fb90759149ecc81a369d15f..35d5706fa21a7fd13d17e19f29c657db3f519538 100644 (file)
--- a/random/random.h
+++ b/random/random.h
@@ -26,7 +26,7 @@ typedef struct s_random {
 #define RAND_STAT_UDEV                 0x04
 #define RAND_STAT_GAUSS                        0x08
 
-#define RAND_BUFSIZE                   (1024) /* 16 mega byte */
+#define RAND_BUFSIZE                   (1024) /* 4 k byte */
 
 #define URAND_MAX                      0xffffffff
 #define URAND_MAX_PLUS_ONE             0x100000000