From 66e48ac1e6b133a607b5fbe2d4d8b7f523659021 Mon Sep 17 00:00:00 2001 From: hackbard Date: Mon, 6 Apr 2009 18:29:01 +0200 Subject: [PATCH] added another diff calc code --- diffusion_calc_ver2.c | 113 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 113 insertions(+) create mode 100644 diffusion_calc_ver2.c diff --git a/diffusion_calc_ver2.c b/diffusion_calc_ver2.c new file mode 100644 index 0000000..ef071a0 --- /dev/null +++ b/diffusion_calc_ver2.c @@ -0,0 +1,113 @@ +/* + * calculation of diffusion coefficient (version 2) + * + * author: frank.zirkelbach@physik.uni-augsburg.de + * + */ + +#define _GNU_SOURCE +#include +#include +#include +#include +#include +#include +#include + +#include "moldyn.h" + +#define CM 1.0e8 + +int usage(char *prog) { + + printf("\nusage:\n"); + printf(" %s \n\n",prog); + + return -1; +} + +int main(int argc,char **argv) { + + t_moldyn moldyn_0,moldyn_1; + int ret,i,a_cnt,b_cnt; + double dx,dy,dz,dc[3]; + + if(argc!=3) { + usage(argv[0]); + return -1; + } + + memset(&moldyn_0,0,sizeof(t_moldyn)); + memset(&moldyn_1,0,sizeof(t_moldyn)); + + a_cnt=0; + b_cnt=0; + + dc[0]=0.0; + dc[1]=0.0; + dc[2]=0.0; + + printf("[diffusion calc] reading save files ...\n"); + ret=moldyn_read_save_file(&moldyn_0,argv[1]); + if(ret) { + printf("[diffusion calc] exit!\n"); + return ret; + } + ret=moldyn_read_save_file(&moldyn_1,argv[2]); + if(ret) { + printf("[diffusion calc] exit!\n"); + return ret; + } + + if(moldyn_0.count!=moldyn_1.count) { + printf("[diffusion calc] atom count mismatch\n"); + return -1; + } + + for(i=0;i0.5) + dx-=0.5; + if(dx<-0.5) + dx+=0.5; + if(dy>0.5) + dy-=0.5; + if(dy<-0.5) + dy+=0.5; + if(dz>0.5) + dz-=0.5; + if(dz<-0.5) + dz+=0.5; + dx*=moldyn_1.dim.x; + dy*=moldyn_1.dim.y; + dz*=moldyn_1.dim.z; + if(moldyn_0.atom[i].brand) { + b_cnt+=1; + dc[1]+=dx*dx+dy*dy+dz*dz; + } + else { + a_cnt+=1; + dc[0]+=dx*dx+dy*dy+dz*dz; + } + dc[2]+=dx*dx+dy*dy+dz*dz; + } + + dc[0]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*a_cnt)); + dc[1]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*b_cnt)); + dc[2]*=(1.0/(6.0*(moldyn_1.time-moldyn_0.time)*moldyn_0.count)); + + dc[0]*=(SECOND/(CM*CM)); + dc[1]*=(SECOND/(CM*CM)); + dc[2]*=(SECOND/(CM*CM)); + + printf("diffusion coefficients: %.10f %.10f %.10f [cm^2/s]\n", + dc[0],dc[1],dc[2]); + + return 0; +} + -- 2.39.2