]> hackdaworld.org Git - physik/posic.git/commitdiff
added fluctuation calc code
authorhackbard <hackbard>
Thu, 5 Jul 2007 14:58:10 +0000 (14:58 +0000)
committerhackbard <hackbard>
Thu, 5 Jul 2007 14:58:10 +0000 (14:58 +0000)
Makefile
fluctuation_calc.c [new file with mode: 0644]
moldyn.c
moldyn.h
random/random.c
sic.c

index bb520c4a81c214d234dbf00d18a3267ca24a0aa5..ce40484dd9feba49a6b9f10ddacd06b293effa79 100644 (file)
--- a/Makefile
+++ b/Makefile
@@ -1,5 +1,4 @@
-CC=gcc-3.4
-#CC=gcc
+CC=gcc
 CFLAGS=-Wall
 CFLAGS+=-O3
 CFLAGS+=-ffloat-store
diff --git a/fluctuation_calc.c b/fluctuation_calc.c
new file mode 100644 (file)
index 0000000..210f382
--- /dev/null
@@ -0,0 +1,121 @@
+/*
+ * calcultae fluctuations and related values
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#define AMU            1.6605388628e-27
+#define METER          1e10
+#define SECOND         1e15
+#define KILOGRAM       (1.0/AMU)
+#define NEWTON         (METER*KILOGRAM/(SECOND*SECOND))
+#define JOULE          (NEWTON*METER)
+#define K_BOLTZMANN    (1.380650524e-23*METER*NEWTON)
+#define K_B2           (K_BOLTZMANN*K_BOLTZMANN)
+#define EV             (1.6021765314e-19*METER*NEWTON)
+
+int get_line(int fd,char *line,int max) {
+
+        int count,ret;
+
+        count=0;
+
+        while(1) {
+                if(count==max) return count;
+                ret=read(fd,line+count,1);
+                if(ret<=0) return ret;
+                if(line[count]=='\n') {
+                        line[count]='\0';
+                        return count+1;
+                }
+                count+=1;
+        }
+}
+
+int main(int argc,char **argv) {
+
+       double from,to;
+       int fd;
+       char file[256];
+       char buf[256];
+       double e,e2,de2,val,time;
+       char *ptr;
+       int i,count,lnr;
+       double np,m,t;
+
+       if(argc<5) {
+               printf("usage:\n");
+               printf("%s <file> <from> <to> <linenumber>\n",argv[0]);
+               printf("\n");
+               printf("optional add: <nr particles> <mass> <temperature>");
+               printf("\n");
+               return -1;
+       }
+       
+       strncpy(file,argv[1],255);
+       from=atof(argv[2]);
+       to=atof(argv[3]);
+       lnr=atoi(argv[4]);
+       np=0;
+       m=0.0;
+       t=0.0;
+
+       if(argc==8) {
+               np=atoi(argv[5]);
+               m=atof(argv[6])*AMU*np;
+               t=atof(argv[7])+273.0;
+       }
+
+       fd=open(file,O_RDONLY);
+       if(fd<2) {
+               perror("fd open");
+               return -1;
+       }
+
+       count=0;
+       e=0.0;
+       e2=0.0;
+
+       while(1) {
+               if(get_line(fd,buf,255)<=0)
+                       break;
+               ptr=strtok(buf," ");
+               time=atof(ptr);
+               if((time>=from)&(time<=to)) {
+                       count+=1;
+                       for(i=1;i<lnr;i++)
+                               ptr=strtok(NULL," ");
+                       val=atof(ptr);
+                       e+=val;
+                       e2+=(val*val);
+               }
+       }
+
+       de2=e2/count-e*e/(count*count);
+       printf("--> fluctuation [eV/atom]: %.20f\n",de2);
+       if(argc==8) {
+               de2*=(np*np)*(EV*EV);
+               printf("    specific heat capacity (T=%f K) [J/(kg K)]:\n",
+                      t);
+               printf("    nve: %f\n",
+                      (1.5*np*K_BOLTZMANN/(1.0-de2/(1.5*np*K_B2*t*t))/
+                      (m*JOULE)));
+               printf("    nvt: %f\n",
+                      (de2/(K_BOLTZMANN*t*t)+1.5*np*K_BOLTZMANN)/
+                      (m*JOULE));
+       }
+
+       close(fd);
+
+       return 0;
+}
+
index cd2a80309c6cf9de56e4b23132debbe9e46b02cd..16c811ed5e0001d905fb7760239fa0b71a055305 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -218,6 +218,14 @@ int set_potential_params(t_moldyn *moldyn,void *params) {
        return 0;
 }
 
+int set_mean_skip(t_moldyn *moldyn,int skip) {
+
+       printf("[moldyn] skip %d steps before starting average calc\n",skip);
+       moldyn->mean_skip=skip;
+
+       return 0;
+}
+
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 
        strncpy(moldyn->vlsdir,dir,127);
@@ -728,8 +736,12 @@ double temperature_calc(t_moldyn *moldyn) {
        /* assume up to date kinetic energy, which is 3/2 N k_B T */
 
        moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+
+       if(moldyn->total_steps<moldyn->mean_skip)
+               return 0;
+
        moldyn->t_sum+=moldyn->t;
-       moldyn->mean_t=moldyn->t_sum/moldyn->total_steps;
+       moldyn->mean_t=moldyn->t_sum/(moldyn->total_steps+1-moldyn->mean_skip);
 
        return moldyn->t;
 }
@@ -826,41 +838,54 @@ double pressure_calc(t_moldyn *moldyn) {
 
        /* virial sum and mean virial */
        moldyn->virial_sum+=v;
-       moldyn->mean_v=moldyn->virial_sum/moldyn->total_steps;
+       if(moldyn->total_steps>=moldyn->mean_skip)
+               moldyn->mean_v=moldyn->virial_sum/
+                              (moldyn->total_steps+1-moldyn->mean_skip);
 
        /* assume up to date kinetic energy */
        moldyn->p=2.0*moldyn->ekin+moldyn->mean_v;
        moldyn->p/=(3.0*moldyn->volume);
-       moldyn->p_sum+=moldyn->p;
-       moldyn->mean_p=moldyn->p_sum/moldyn->total_steps;
+       if(moldyn->total_steps>=moldyn->mean_skip) {
+               moldyn->p_sum+=moldyn->p;
+               moldyn->mean_p=moldyn->p_sum/
+                              (moldyn->total_steps+1-moldyn->mean_skip);
+       }
 
        /* pressure from 'absolute coordinates' virial */
        virial=&(moldyn->virial);
        v=virial->xx+virial->yy+virial->zz;
        moldyn->gp=2.0*moldyn->ekin+v;
        moldyn->gp/=(3.0*moldyn->volume);
-       moldyn->gp_sum+=moldyn->gp;
-       moldyn->mean_gp=moldyn->gp_sum/moldyn->total_steps;
+       if(moldyn->total_steps>=moldyn->mean_skip) {
+               moldyn->gp_sum+=moldyn->gp;
+               moldyn->mean_gp=moldyn->gp_sum/
+                              (moldyn->total_steps+1-moldyn->mean_skip);
+       }
 
        return moldyn->p;
 }
 
 int energy_fluctuation_calc(t_moldyn *moldyn) {
 
+       if(moldyn->total_steps<moldyn->mean_skip)
+               return 0;
+
        /* assume up to date energies */
 
        /* kinetic energy */
        moldyn->k_sum+=moldyn->ekin;
        moldyn->k2_sum+=(moldyn->ekin*moldyn->ekin);
-       moldyn->k_mean=moldyn->k_sum/moldyn->total_steps;
-       moldyn->k2_mean=moldyn->k2_sum/moldyn->total_steps;
+       moldyn->k_mean=moldyn->k_sum/(moldyn->total_steps+1-moldyn->mean_skip);
+       moldyn->k2_mean=moldyn->k2_sum/
+                       (moldyn->total_steps+1-moldyn->mean_skip);
        moldyn->dk2_mean=moldyn->k2_mean-(moldyn->k_mean*moldyn->k_mean);
 
        /* potential energy */
        moldyn->v_sum+=moldyn->energy;
        moldyn->v2_sum+=(moldyn->energy*moldyn->energy);
-       moldyn->v_mean=moldyn->v_sum/moldyn->total_steps;
-       moldyn->v2_mean=moldyn->v2_sum/moldyn->total_steps;
+       moldyn->v_mean=moldyn->v_sum/(moldyn->total_steps+1-moldyn->mean_skip);
+       moldyn->v2_mean=moldyn->v2_sum/
+                       (moldyn->total_steps+1-moldyn->mean_skip);
        moldyn->dv2_mean=moldyn->v2_mean-(moldyn->v_mean*moldyn->v_mean);
 
        return 0;
@@ -870,6 +895,10 @@ int get_heat_capacity(t_moldyn *moldyn) {
 
        double temp2,ighc;
 
+       /* averages needed for heat capacity calc */
+       if(moldyn->total_steps<moldyn->mean_skip)
+               return 0;
+
        /* (temperature average)^2 */
        temp2=moldyn->mean_t*moldyn->mean_t;
        printf("[moldyn] specific heat capacity for T=%f K [J/(kg K)]\n",
@@ -890,6 +919,7 @@ int get_heat_capacity(t_moldyn *moldyn) {
 
        printf("  NVE: %f\n",moldyn->c_v_nve*KILOGRAM/JOULE);
        printf("  NVT: %f\n",moldyn->c_v_nvt*KILOGRAM/JOULE);
+printf("  --> <dV2> sim: %f experimental: %f\n",moldyn->dv2_mean,1.5*moldyn->count*K_B2*moldyn->mean_t*moldyn->mean_t*(1.0-1.5*moldyn->count*K_BOLTZMANN/(700*moldyn->mass*JOULE/KILOGRAM)));
 
        return 0;
 }
index 2d1a438b3935b05cc8dfd791a6505b8abf03a093..868db6267784faa8d684d265e6a976668776f10f 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -102,6 +102,8 @@ typedef struct s_moldyn {
 
        t_linkcell lc;          /* linked cell list interface */
 
+       int mean_skip;          /* amount of steps without average calc */
+
        double t_ref;           /* reference temperature */
        double t;               /* actual temperature */
        double t_sum;           /* sum over all t */
@@ -309,7 +311,7 @@ typedef struct s_moldyn {
 
 #define TM_CHI_SIC             0.9776
 
-#define TM_LC_3C_SIC           (0.432e-9*METER)                /* A */
+#define TM_LC_SIC              4.32                            /* A */
 
 #define ALBE_R_SI              (2.82-0.14)
 #define ALBE_S_SI              (2.82+0.14)
@@ -323,7 +325,7 @@ typedef struct s_moldyn {
 #define ALBE_D_SI              0.81472
 #define ALBE_H_SI              0.259
 
-#define LC_SI_ALBE             5.429
+#define ALBE_LC_SI             5.429
 
 #define ALBE_R_C               (2.00-0.15)
 #define ALBE_S_C               (2.00+0.15)
@@ -337,7 +339,7 @@ typedef struct s_moldyn {
 #define ALBE_D_C               6.28433
 #define ALBE_H_C               0.5556
 
-#define LC_C_ALBE              3.566
+#define ALBE_LC_C              3.566
 
 #define ALBE_R_SIC             (2.40-0.20)
 #define ALBE_S_SIC             (2.40+0.10)
@@ -351,7 +353,7 @@ typedef struct s_moldyn {
 #define ALBE_D_SIC             180.314
 #define ALBE_H_SIC             0.68
 
-#define LC_SIC_ALBE            4.359
+#define ALBE_LC_SIC            4.359
 
 
 /*
@@ -393,6 +395,8 @@ int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func);
 int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func);
 int set_potential_params(t_moldyn *moldyn,void *params);
 
+int set_mean_skip(t_moldyn *moldyn,int skip);
+
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
 int moldyn_set_report(t_moldyn *moldyn,char *author,char *title);
 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
index df779809b5a34b627eb6ccfe35ad03835f709952..902cdef4738d9d95f81b58b5fc72333948a60a1b 100644 (file)
@@ -36,6 +36,7 @@ int rand_init(t_random *random,char *randomfile,int logfd) {
                        return -1;
                }
        }
+       random->buffer=NULL;
        random->buffer=malloc(RAND_BUFSIZE*sizeof(unsigned int));
        if(random->buffer==NULL) {
                perror("malloc random buffer");
diff --git a/sic.c b/sic.c
index 426998a15fefbf87e436408189a3a1a746f1d8eb..30de32874025ff60b8dfae01224730c95032c037 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -210,12 +210,13 @@ int main(int argc,char **argv) {
 
        /* set (initial) dimensions of simulation volume */
 #ifdef ALBLE
-       set_dim(&md,6*LC_SI_ALBE,6*LC_SI_ALBE,6*LC_SI_ALBE,TRUE);
-       //set_dim(&md,6*LC_C_ALBE,6*LC_C_ALBE,6*LC_C_ALBE,TRUE);
-       //set_dim(&md,6*LC_SIC_ALBE,6*LC_SIC_ALBE,6*LC_SIC_ALBE,TRUE);
+       //set_dim(&md,8*ALBE_LC_SI,8*ALBE_LC_SI,8*ALBE_LC_SI,TRUE);
+       //set_dim(&md,8*ALBE_LC_C,8*ALBE_LC_C,8*ALBE_LC_C,TRUE);
+       set_dim(&md,8*ALBE_LC_SIC,8*ALBE_LC_SIC,8*ALBE_LC_SIC,TRUE);
 #else
-       set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
-       //set_dim(&md,6*LC_C,6*LC_C,6*LC_C,TRUE);
+       //set_dim(&md,8*LC_SI,8*LC_SI,8*LC_SI,TRUE);
+       //set_dim(&md,8*LC_C,8*LC_C,8*LC_C,TRUE);
+       set_dim(&md,8*TM_LC_SIC,8*TM_LC_SIC,8*TM_LC_SIC,TRUE);
 #endif
 
        /* set periodic boundary conditions in all directions */
@@ -223,28 +224,40 @@ int main(int argc,char **argv) {
 
        /* create the lattice / place atoms */
 #ifdef ALBLE
-       create_lattice(&md,DIAMOND,LC_SI_ALBE,SI,M_SI,
+       //create_lattice(&md,DIAMOND,ALBE_LC_SI,SI,M_SI,
+       //create_lattice(&md,DIAMOND,ALBE_LC_C,C,M_C,
 #else
-       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
-       //create_lattice(&md,DIAMOND,LC_C_ALBE,C,M_C,
+       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
 #endif
-                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+       //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
        //               ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,6,6,6,NULL);
-       //               1,6,6,6,NULL);
+       //               0,8,8,8,NULL);
+       //               1,8,8,8,NULL);
 
-       /* create centered zinc blende lattice */
-       /*
-       r.x=0.5*0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x;
-       create_lattice(&md,FCC,LC_SIC_ALBE,SI,M_SI,
+       /* create zinkblende structure */
+       /**/
+#ifdef ALBE
+       r.x=0.5*0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,ALBE_LC_SIC,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      0,8,8,8,&r);
+       r.x+=0.25*ALBE_LC_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,ALBE_LC_SIC,C,M_C,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-                      0,6,6,6,&r);
-       r.x+=0.25*LC_SIC_ALBE; r.y=r.x; r.z=r.x;
-       create_lattice(&md,FCC,LC_SIC_ALBE,C,M_C,
+                      1,8,8,8,&r);
+#else
+       r.x=0.5*0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,TM_LC_SIC,SI,M_SI,
+                      ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                      0,8,8,8,&r);
+       r.x+=0.25*TM_LC_SIC; r.y=r.x; r.z=r.x;
+       create_lattice(&md,FCC,TM_LC_SIC,C,M_C,
                       ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
-                      1,6,6,6,&r);
-       */
+                      1,8,8,8,&r);
+#endif
+       /**/
 
+       /* check for right atom placing */
        moldyn_bc_check(&md);
 
        /* testing configuration */
@@ -283,6 +296,9 @@ int main(int argc,char **argv) {
        set_temperature(&md,atof(argv[2])+273.0);
        set_pressure(&md,BAR);
 
+       /* set amount of steps to skip before average calc */
+       set_mean_skip(&md,500);
+
        /* set p/t scaling */
        //set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
        //set_pt_scale(&md,P_SCALE_BERENDSEN,0.001,
@@ -295,7 +311,7 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        /* initial configuration */
-       moldyn_add_schedule(&md,3000,1.0);
+       moldyn_add_schedule(&md,1500,1.0);
        //moldyn_add_schedule(&md,1000,1.0);
        //moldyn_add_schedule(&md,1000,1.0);
        //moldyn_add_schedule(&md,1000,1.0);