diffusion_calc_ver2.c

   1 /*
   2  * calculation of diffusion coefficient (version 2)
   3  *
   4  * author: frank.zirkelbach@physik.uni-augsburg.de
   5  *
   6  */
   7
   8 #define _GNU_SOURCE
   9 #include <stdio.h>
  10 #include <stdlib.h>
  11 #include <unistd.h>
  12 #include <string.h>
  13 #include <sys/types.h>
  14 #include <sys/stat.h>
  15 #include <fcntl.h>
  16
  17 #include "moldyn.h"
  18
  19 #define CM      1.0e8
  20
  21 int usage(char *prog) {
  22
  23         printf("\nusage:\n");
  24         printf("  %s <save file 1> <save file 2>\n\n",prog);
  25
  26         return -1;
  27 }
  28
  29 int main(int argc,char **argv) {
  30
  31         t_moldyn moldyn_0,moldyn_1;
  32         int ret,i,a_cnt,b_cnt;
  33         double dx,dy,dz,dc[3];
  34
  35         if(argc!=3) {
  36                 usage(argv[0]);
  37                 return -1;
  38         }
  39
  40         memset(&moldyn_0,0,sizeof(t_moldyn));
  41         memset(&moldyn_1,0,sizeof(t_moldyn));
  42
  43         a_cnt=0;
  44         b_cnt=0;
  45
  46         dc[0]=0.0;
  47         dc[1]=0.0;
  48         dc[2]=0.0;
  49
  50         printf("[diffusion calc] reading save files ...\n");
  51         ret=moldyn_read_save_file(&moldyn_0,argv[1]);
  52         if(ret) {
  53                 printf("[diffusion calc] exit!\n");
  54                 return ret;
  55         }
  56         ret=moldyn_read_save_file(&moldyn_1,argv[2]);
  57         if(ret) {
  58                 printf("[diffusion calc] exit!\n");
  59                 return ret;
  60         }
  61
  62         if(moldyn_0.count!=moldyn_1.count) {
  63                 printf("[diffusion calc] atom count mismatch\n");
  64                 return -1;
  65         }
  66
  67         for(i=0;i<moldyn_0.count;i++) {
  68                 dx=moldyn_0.atom[i].r.x/moldyn_0.dim.x;
  69                 dy=moldyn_0.atom[i].r.y/moldyn_0.dim.y;
  70                 dz=moldyn_0.atom[i].r.z/moldyn_0.dim.z;
  71                 dx-=(moldyn_1.atom[i].r.x/moldyn_1.dim.x);
  72                 dy-=(moldyn_1.atom[i].r.y/moldyn_1.dim.y);
  73                 dz-=(moldyn_1.atom[i].r.z/moldyn_1.dim.z);
  74                 if(dx>0.5)
  75                         dx-=0.5;
  76                 if(dx<-0.5)
  77                         dx+=0.5;
  78                 if(dy>0.5)
  79                         dy-=0.5;
  80                 if(dy<-0.5)
  81                         dy+=0.5;
  82                 if(dz>0.5)
  83                         dz-=0.5;
  84                 if(dz<-0.5)
  85                         dz+=0.5;
  86                 dx*=moldyn_1.dim.x;
  87                 dy*=moldyn_1.dim.y;
  88                 dz*=moldyn_1.dim.z;
  89                 if(moldyn_0.atom[i].brand) {
  90                         b_cnt+=1;
  91                         dc[1]+=dx*dx+dy*dy+dz*dz;
  92                 }
  93                 else {
  94                         a_cnt+=1;
  95                         dc[0]+=dx*dx+dy*dy+dz*dz;
  96                 }
  97                 dc[2]+=dx*dx+dy*dy+dz*dz;
  98         }
  99
 100         dc[0]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*a_cnt));
 101         dc[1]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*b_cnt));
 102         dc[2]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*moldyn_0.count));
 103
 104         dc[0]*=(SECOND/(CM*CM));
 105         dc[1]*=(SECOND/(CM*CM));
 106         dc[2]*=(SECOND/(CM*CM));
 107
 108         printf("diffusion coefficients: %.10f %.10f %.10f [cm^2/s]\n",
 109                dc[0],dc[1],dc[2]);
 110
 111         return 0;
 112 }
 113