]> hackdaworld.org Git - physik/posic.git/blob - fluctuation_calc.c
optional elment and brand type in fill cmd now
[physik/posic.git] / 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