From 2346d076f3d0955fb3a7bfb47f2b17b9eae6dd5b Mon Sep 17 00:00:00 2001 From: hackbard Date: Wed, 9 Apr 2008 14:47:25 +0200 Subject: [PATCH] introduce a 2body bond function with callback, modified pair corr calc approp + introduced minimal bond analyze code --- Makefile | 9 +- bond_analyze.c | 54 ++++++ moldyn.c | 357 +++++++++++++++++++++++++++------------- moldyn.h | 22 +++ pair_corr_calc_script | 6 +- pair_correlation_calc.c | 2 + 6 files changed, 329 insertions(+), 121 deletions(-) create mode 100644 bond_analyze.c diff --git a/Makefile b/Makefile index 1b875fe..c839759 100644 --- a/Makefile +++ b/Makefile @@ -22,7 +22,10 @@ DEPS = moldyn.o random/random.o list/list.o DEPS += potentials/lennard_jones.o potentials/harmonic_oscillator.o DEPS += potentials/tersoff.o potentials/albe.o -all: posic sic fluctuation_calc postproc pair_correlation_calc diffusion_calc +ALL = posic sic fluctuation_calc postproc pair_correlation_calc diffusion_calc +ALL += bond_analyze + +all: $(ALL) posic: $(DEPS) @@ -34,6 +37,8 @@ pair_correlation_calc: $(DEPS) diffusion_calc: $(DEPS) +bond_analyze: $(DEPS) + .PHONY:clean clean: - rm -vf sic postproc fluctuation_calc pair_correlation_calc *.o */*.o + rm -vf $(ALL) *.o */*.o diff --git a/bond_analyze.c b/bond_analyze.c new file mode 100644 index 0000000..0feffef --- /dev/null +++ b/bond_analyze.c @@ -0,0 +1,54 @@ +/* + * bonding analyzation code + * + * author: frank.zirkelbach@physik.uni-augsburg.de + * + */ + +#define _GNU_SOURCE +#include +//#include +//#include +//#include +//#include +//#include +//#include + +#include "moldyn.h" + +int usage(char *prog) { + + printf("\nusage:\n"); + printf(" %s \n\n",prog); + + return -1; +} + +int main(int argc,char **argv) { + + t_moldyn moldyn; + int ret; + double quality; + + if(argc!=2) { + usage(argv[0]); + return -1; + } + + memset(&moldyn,0,sizeof(t_moldyn)); + + printf("[bond analyze] reading save file ...\n"); + ret=moldyn_read_save_file(&moldyn,argv[1]); + if(ret) { + printf("[bond analyze] exit!\n"); + return ret; + } + + bond_analyze(&moldyn,&quality); + + printf("[bond analyze] quality = %f\n",quality); + + moldyn_free_save_file(&moldyn); + + return 0; +} diff --git a/moldyn.c b/moldyn.c index 91fead7..84cf700 100644 --- a/moldyn.c +++ b/moldyn.c @@ -2323,6 +2323,13 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { return 0; } +int moldyn_free_save_file(t_moldyn *moldyn) { + + free(moldyn->atom); + + return 0; +} + int moldyn_load(t_moldyn *moldyn) { // later ... @@ -2330,6 +2337,78 @@ int moldyn_load(t_moldyn *moldyn) { return 0; } +/* + * function to find/callback all combinations of 2 body bonds + */ + +int process_2b_bonds(t_moldyn *moldyn,void *data, + int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc)) { + + t_linkcell *lc; +#ifdef STATIC_LISTS + int *neighbour[27]; + int p; +#else + t_list neighbour[27]; +#endif + u8 bc; + t_atom *itom,*jtom; + int i,j; + t_list *this; + + lc=&(moldyn->lc); + + link_cell_init(moldyn,VERBOSE); + + itom=moldyn->atom; + + for(i=0;icount;i++) { + /* neighbour indexing */ + link_cell_neighbour_index(moldyn, + (itom[i].r.x+moldyn->dim.x/2)/lc->x, + (itom[i].r.y+moldyn->dim.y/2)/lc->x, + (itom[i].r.z+moldyn->dim.z/2)/lc->x, + neighbour); + + for(j=0;j<27;j++) { + + bc=(jdnlc)?0:1; + +#ifdef STATIC_LISTS + p=0; + + while(neighbour[j][p]!=0) { + + jtom=&(moldyn->atom[neighbour[j][p]]); + p++; +#else + this=&(neighbour[j]); + list_reset_f(this); + + if(this->start==NULL) + continue; + + do { + + jtom=this->current->data; +#endif + + /* process bond */ + process(moldyn,&(itom[i]),jtom,data,bc); + +#ifdef STATIC_LISTS + } +#else + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); +#endif + } + } + + return 0; + +} + /* * post processing functions */ @@ -2399,33 +2478,80 @@ int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) { return 0; } -int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { +int bonding_analyze(t_moldyn *moldyn,double *cnt) { + + return 0; +} + +int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, + t_atom *jtom,void *data,u8 bc) { - int slots; - double *stat; - int i,j; - t_linkcell *lc; -#ifdef STATIC_LISTS - int *neighbour[27]; - int p; -#else - t_list neighbour[27]; -#endif - t_atom *itom,*jtom; - t_list *this; - unsigned char bc; t_3dvec dist; double d; - double norm; - int o,s; - unsigned char ibrand; + int s; + t_pcc *pcc; - lc=&(moldyn->lc); + /* only count pairs once, + * skip same atoms */ + if(itom->tag>=jtom->tag) + return 0; + + /* + * pair correlation calc + */ + + /* get pcc data */ + pcc=data; + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + + /* ignore if greater or equal 2 times cutoff */ + if(d>=4.0*moldyn->cutoff_square) + return 0; + + /* fill the slots */ + d=sqrt(d); + s=(int)(d/pcc->dr); + + /* should never happen but it does 8) - + * related to -ffloat-store problem! */ + if(s>=pcc->o1) { + printf("[moldyn] WARNING: pcc (%d/%d)", + s,pcc->o1); + printf("\n"); + s=pcc->o1-1; + } + + if(itom->brand!=jtom->brand) { + /* mixed */ + pcc->stat[s]+=1; + } + else { + /* type a - type a bonds */ + if(itom->brand==0) + pcc->stat[s+pcc->o1]+=1; + else + /* type b - type b bonds */ + pcc->stat[s+pcc->o2]+=1; + } - slots=2.0*moldyn->cutoff/dr; - o=2*slots; + return 0; +} - if(slots*dr<=moldyn->cutoff) +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { + + t_pcc pcc; + double norm; + int i; + + pcc.dr=dr; + pcc.o1=2.0*moldyn->cutoff/dr; + pcc.o2=2*pcc.o1; + + if(pcc.o1*dr<=moldyn->cutoff) printf("[moldyn] WARNING: pcc (low #slots)\n"); printf("[moldyn] pair correlation calc info:\n"); @@ -2435,134 +2561,129 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg); if(ptr!=NULL) { - stat=(double *)ptr; + pcc.stat=(double *)ptr; } else { - stat=(double *)malloc(3*slots*sizeof(double)); - if(stat==NULL) { + pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double)); + if(pcc.stat==NULL) { perror("[moldyn] pair correlation malloc"); return -1; } } - memset(stat,0,3*slots*sizeof(double)); + memset(pcc.stat,0,3*pcc.o1*sizeof(double)); - link_cell_init(moldyn,VERBOSE); + /* process */ + process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process); - itom=moldyn->atom; - - for(i=0;icount;i++) { - /* neighbour indexing */ - link_cell_neighbour_index(moldyn, - (itom[i].r.x+moldyn->dim.x/2)/lc->x, - (itom[i].r.y+moldyn->dim.y/2)/lc->x, - (itom[i].r.z+moldyn->dim.z/2)/lc->x, - neighbour); + /* normalization */ + for(i=1;i 2 pi r r dr + // ... and actually it's a constant times r^2 + norm=i*i*dr*dr; + pcc.stat[i]/=norm; + pcc.stat[pcc.o1+i]/=norm; + pcc.stat[pcc.o2+i]/=norm; + } + /* */ - /* brand of atom i */ - ibrand=itom[i].brand; - - for(j=0;j<27;j++) { + if(ptr==NULL) { + /* todo: store/print pair correlation function */ + free(pcc.stat); + } - bc=(jdnlc)?0:1; + return 0; +} -#ifdef STATIC_LISTS - p=0; +int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc) { - while(neighbour[j][p]!=0) { + t_ba *ba; + t_3dvec dist; + double d; - jtom=&(moldyn->atom[neighbour[j][p]]); - p++; -#else - this=&(neighbour[j]); - list_reset_f(this); + if(itom->tag>=jtom->tag) + return 0; - if(this->start==NULL) - continue; + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom->r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); - do { + /* ignore if greater or equal cutoff */ + if(d>=moldyn->cutoff_square) + return 0; - jtom=this->current->data; -#endif - /* only count pairs once, - * skip same atoms */ - if(itom[i].tag>=jtom->tag) - continue; + /* now count this bonding ... */ + ba=data; - /* - * pair correlation calc - */ + /* increase total bond counter + * ... double counting! + */ + ba->tcnt+=2; - /* distance */ - v3_sub(&dist,&(jtom->r),&(itom[i].r)); - if(bc) check_per_bound(moldyn,&dist); - d=v3_absolute_square(&dist); + if(itom->brand==0) + ba->acnt[jtom->tag]+=1; + else + ba->bcnt[jtom->tag]+=1; + + if(jtom->brand==0) + ba->acnt[itom->tag]+=1; + else + ba->bcnt[itom->tag]+=1; - /* ignore if greater or equal cutoff */ - if(d>=4.0*moldyn->cutoff_square) - continue; + return 0; +} - /* fill the slots */ - d=sqrt(d); - s=(int)(d/dr); - - /* should never happen but it does 8) - - * related to -ffloat-store problem! */ - if(s>=slots) { - printf("[moldyn] WARNING: pcc (%d/%d)", - s,slots); - printf("\n"); - s=slots-1; - } +int bond_analyze(t_moldyn *moldyn,double *quality) { - if(ibrand!=jtom->brand) { - /* mixed */ - stat[s]+=1; - } - else { - /* type a - type a bonds */ - if(ibrand==0) - stat[s+slots]+=1; - else - /* type b - type b bonds */ - stat[s+o]+=1; - } -#ifdef STATIC_LISTS - } -#else - } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); -#endif - } - } + // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total - /* normalization */ - for(i=1;i 2 pi r r dr - // ... and actually it's a constant times r^2 - norm=i*i*dr*dr; - stat[i]/=norm; - stat[slots+i]/=norm; - stat[o+i]/=norm; + int qcnt; + t_ba ba; + int i; + t_atom *atom; + + ba.acnt=malloc(moldyn->count*sizeof(int)); + if(ba.acnt==NULL) { + perror("[moldyn] bond analyze malloc (a)"); + return -1; } - /* */ + memset(ba.acnt,0,moldyn->count*sizeof(int)); - if(ptr==NULL) { - /* todo: store/print pair correlation function */ - free(stat); + ba.bcnt=malloc(moldyn->count*sizeof(int)); + if(ba.bcnt==NULL) { + perror("[moldyn] bond analyze malloc (b)"); + return -1; } + memset(ba.bcnt,0,moldyn->count*sizeof(int)); - free(moldyn->atom); + ba.tcnt=0; - link_cell_shutdown(moldyn); + qcnt=0; - return 0; -} + atom=moldyn->atom; -int analyze_bonds(t_moldyn *moldyn) { + process_2b_bonds(moldyn,&ba,bond_analyze_process); - - + for(i=0;icount;i++) { + if(atom[i].brand==0) { + if((ba.acnt[i]==0)&(ba.bcnt[i]==4)) + qcnt+=4; + } + else { + if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) + qcnt+=4; + } + } + +printf("%d %d\n",qcnt,ba.tcnt); + if(quality) + *quality=1.0*qcnt/ba.tcnt; + else + printf("[moldyn] bond analyze: quality = %f\n", + 1.0*qcnt/ba.tcnt); return 0; } diff --git a/moldyn.h b/moldyn.h index 83b3a86..1647c1c 100644 --- a/moldyn.h +++ b/moldyn.h @@ -219,6 +219,19 @@ typedef struct s_moldyn { double debug; /* debugging stuff, ignore */ } t_moldyn; +typedef struct s_pcc { + int o1; + int o2; + double dr; + double *stat; +} t_pcc; + +typedef struct s_ba { + int *acnt; + int *bcnt; + int tcnt; +} t_ba; + /* * * defines @@ -494,12 +507,21 @@ int check_per_bound(t_moldyn *moldyn,t_3dvec *a); int moldyn_bc_check(t_moldyn *moldyn); int moldyn_read_save_file(t_moldyn *moldyn,char *file); +int moldyn_free_save_file(t_moldyn *moldyn); int moldyn_load(t_moldyn *moldyn); +int process_2b_bonds(t_moldyn *moldyn,void *data, + int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc)); int get_line(int fd,char *line,int max); int pair_correlation_init(t_moldyn *moldyn,double dr); int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc); +int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, + t_atom *jtom,void *data,u8 bc); int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr); +int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc); +int bond_analyze(t_moldyn *moldyn,double *quality); int visual_init(t_moldyn *moldyn,char *filebase); int visual_atoms(t_moldyn *moldyn); diff --git a/pair_corr_calc_script b/pair_corr_calc_script index 7f25a74..b005960 100755 --- a/pair_corr_calc_script +++ b/pair_corr_calc_script @@ -27,7 +27,11 @@ fi if [ "$3" = "g" ]; then -pdir=`dirname $1` +if [ -d $1 ]; then + pdir=$1 +else + pdir=`dirname $1` +fi pfile=$pdir/pair_corr.scr cat > $pfile <<-EOF diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c index e5903aa..c7716a3 100644 --- a/pair_correlation_calc.c +++ b/pair_correlation_calc.c @@ -76,5 +76,7 @@ int main(int argc,char **argv) { free(stat); + moldyn_free_save_file(&moldyn); + return 0; } -- 2.39.2