fluctuation_calc.c

   1 /*
   2  * calcultae fluctuations and related values
   3  *
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10 #include <unistd.h>
  11 #include <string.h>
  12 #include <sys/types.h>
  13 #include <sys/stat.h>
  14 #include <fcntl.h>
  15
  16 #define AMU             1.6605388628e-27
  17 #define METER           1e10
  18 #define SECOND          1e15
  19 #define KILOGRAM        (1.0/AMU)
  20 #define NEWTON          (METER*KILOGRAM/(SECOND*SECOND))
  21 #define JOULE           (NEWTON*METER)
  22 #define K_BOLTZMANN     (1.380650524e-23*METER*NEWTON)
  23 #define K_B2            (K_BOLTZMANN*K_BOLTZMANN)
  24 #define EV              (1.6021765314e-19*METER*NEWTON)
  25
  26 int get_line(int fd,char *line,int max) {
  27
  28         int count,ret;
  29
  30         count=0;
  31
  32         while(1) {
  33                 if(count==max) return count;
  34                 ret=read(fd,line+count,1);
  35                 if(ret<=0) return ret;
  36                 if(line[count]=='\n') {
  37                         line[count]='\0';
  38                         return count+1;
  39                 }
  40                 count+=1;
  41         }
  42 }
  43
  44 int main(int argc,char **argv) {
  45
  46         double from,to;
  47         int fd;
  48         char file[256];
  49         char buf[256];
  50         double e,e2,de2,val,time;
  51         char *ptr;
  52         int i,count,lnr;
  53         double np,m,t;
  54
  55         if(argc<5) {
  56                 printf("usage:\n");
  57                 printf("%s <file> <from> <to> <linenumber>\n",argv[0]);
  58                 printf("\n");
  59                 printf("optional add: <nr particles> <mass> <temperature>");
  60                 printf("\n");
  61                 return -1;
  62         }
  63
  64         strncpy(file,argv[1],255);
  65         from=atof(argv[2]);
  66         to=atof(argv[3]);
  67         lnr=atoi(argv[4]);
  68         np=0;
  69         m=0.0;
  70         t=0.0;
  71
  72         if(argc==8) {
  73                 np=atoi(argv[5]);
  74                 m=atof(argv[6])*AMU*np;
  75                 t=atof(argv[7])+273.0;
  76         }
  77
  78         fd=open(file,O_RDONLY);
  79         if(fd<2) {
  80                 perror("fd open");
  81                 return -1;
  82         }
  83
  84         count=0;
  85         e=0.0;
  86         e2=0.0;
  87
  88         while(1) {
  89                 if(get_line(fd,buf,255)<=0)
  90                         break;
  91                 ptr=strtok(buf," ");
  92                 time=atof(ptr);
  93                 if((time>=from)&(time<=to)) {
  94                         count+=1;
  95                         for(i=1;i<lnr;i++)
  96                                 ptr=strtok(NULL," ");
  97                         val=atof(ptr);
  98                         e+=val;
  99                         e2+=(val*val);
 100                 }
 101         }
 102
 103         de2=e2/count-e*e/(count*count);
 104         printf("--> fluctuation [eV/atom]: %.20f\n",de2);
 105         if(argc==8) {
 106                 de2*=((np*np*EV*EV));
 107                 printf("    specific heat capacity (T=%f K) [J/(kg K)]:\n",
 108                        t);
 109                 printf("    nve: %f\n",
 110                        (1.5*np*K_BOLTZMANN/(1.0-de2/(1.5*np*K_B2*t*t))/
 111                        (m*JOULE)));
 112                 printf("    nvt: %f\n",
 113                        (de2/(K_BOLTZMANN*t*t)+1.5*np*K_BOLTZMANN)/
 114                        (m*JOULE));
 115         }
 116
 117         close(fd);
 118
 119         return 0;
 120 }
 121