From: hackbard Date: Thu, 24 Jan 2008 16:39:00 +0000 (+0100) Subject: added pair correlation calc code X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=commitdiff_plain;h=dff2917d22ed07707d222bc10fab7370356699dc added pair correlation calc code --- diff --git a/Makefile b/Makefile index 0e7b122..73645df 100644 --- a/Makefile +++ b/Makefile @@ -17,7 +17,7 @@ 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 +all: posic sic fluctuation_calc postproc pair_correlation_calc posic: $(DEPS) @@ -25,6 +25,8 @@ sic: $(DEPS) postproc: $(DEPS) +pair_correlation_calc: $(DEPS) + .PHONY:clean clean: rm -vf sic postproc fluctuation_calc *.o */*.o diff --git a/moldyn.c b/moldyn.c index ad4f900..2e0e468 100644 --- a/moldyn.c +++ b/moldyn.c @@ -2057,6 +2057,13 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) { } size=moldyn->count*sizeof(t_atom); + + moldyn->atom=(t_atom *)malloc(size); + if(moldyn->atom==NULL) { + perror("[moldyn] load save file malloc (atoms)"); + return -1; + } + cnt=read(fd,moldyn->atom,size); if(cnt!=size) { perror("[moldyn] load save file read (atoms)"); @@ -2103,49 +2110,130 @@ int pair_correlation_init(t_moldyn *moldyn,double dr) { return 0; } -int calculate_pair_correlation(t_moldyn *moldyn,double dr) { +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) { int slots; - int *stat; - int i; + double *stat; + int i,j; t_linkcell *lc; t_list neighbour[27]; t_atom *itom,*jtom; t_list *this; + unsigned char bc; + t_3dvec dist; + double d,norm; + int o,s; + unsigned char ibrand; lc=&(moldyn->lc); slots=(int)(moldyn->cutoff/dr); + o=2*slots; - stat=(int *)malloc(3*slots*sizeof(int)); - if(stat==NULL) { - perror("[moldyn] pair correlation malloc"); - return -1; + if(ptr!=NULL) { + stat=(double *)ptr; + } + else { + stat=(double *)malloc(3*slots*sizeof(double)); + if(stat==NULL) { + perror("[moldyn] pair correlation malloc"); + return -1; + } } + memset(stat,0,3*slots*sizeof(double)); + link_cell_init(moldyn,VERBOSE); - - for(i=0;icells;i++) { - /* check for atoms */ - itom=lc->subcell[i].current->data; - if(itom==NULL) - continue; - /* pick first atom and do neighbour indexing */ + itom=moldyn->atom; + + for(i=0;icount;i++) { + /* neighbour indexing */ link_cell_neighbour_index(moldyn, - (itom->r.x+moldyn->dim.x/2)/lc->x, - (itom->r.y+moldyn->dim.y/2)/lc->x, - (itom->r.z+moldyn->dim.z/2)/lc->x, + (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); - /* calculation of g(r) */ - do { - itom=lc->subcell[i].current->data; -// GO ON HERE ... - } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);); + /* brand of atom i */ + ibrand=itom[i].brand; + + for(j=0;j<27;j++) { + /* prepare the neighbour cell list */ + this=&(neighbour[j]); + list_reset_f(this); + + /* check for atoms */ + if(this->start==NULL) + continue; + + /* boundary check */ + bc=(jdnlc)?0:1; + + do { + jtom=this->current->data; + + + if(jtom==&(itom[i])) + continue; + + /* only count pairs once */ + if(itom[i].tag>jtom->tag) + continue; + + /* + * pair correlation calc + */ + + /* distance */ + v3_sub(&dist,&(jtom->r),&(itom[i].r)); + if(bc) check_per_bound(moldyn,&dist); + d=v3_absolute_square(&dist); + /* ignore if greater cutoff */ + if(d>moldyn->cutoff_square) + continue; + + /* fill the slots */ + d=sqrt(d); + s=(int)(d/dr); + + 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; + } + + } while(list_next_f(this)!=L_NO_NEXT_ELEMENT); + } } + /* normalization */ + for(i=1;i 2 pi r r dr + */ + norm=2*M_PI*moldyn->count*(i*dr*i*dr)*dr; + stat[i]/=norm; + stat[slots+i]/=norm; + stat[o+i]/=norm; + } + + if(ptr==NULL) { + /* todo: store/print pair correlation function */ + free(stat); + } + + free(moldyn->atom); + link_cell_shutdown(moldyn); return 0; diff --git a/moldyn.h b/moldyn.h index 2d4ca48..a7a81e9 100644 --- a/moldyn.h +++ b/moldyn.h @@ -475,8 +475,13 @@ 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_load(t_moldyn *moldyn); int get_line(int fd,char *line,int max); +int pair_correlation_init(t_moldyn *moldyn,double dr); +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr); + 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 new file mode 100755 index 0000000..27dc0ab --- /dev/null +++ b/pair_corr_calc_script @@ -0,0 +1,29 @@ +#!/bin/sh + +# +# calculate pair correlation functions for a whole simulation output +# frank.zirkelbach@physik.uni-augsburg.de +# + +do_it() { + echo "processing $1 ..." + ./pair_correlation_calc $1 $2 + trgab=`echo $1 | sed 's%s-%pair_corr_ab-%' | sed 's%.save%%'` + trgaa=`echo $1 | sed 's%s-%pair_corr_aa-%' | sed 's%.save%%'` + trgbb=`echo $1 | sed 's%s-%pair_corr_bb-%' | sed 's%.save%%'` + mv pair_corr_func_ab.txt $trgab + mv pair_corr_func_aa.txt $trgaa + mv pair_corr_func_bb.txt $trgbb + echo "done" +} + +if [ -d $1 ]; then + for file in $1/*.save; do + do_it $file $2 + done +fi + +if [ -f $1 ]; then + do_it $1 $2 +fi + diff --git a/pair_correlation_calc.c b/pair_correlation_calc.c new file mode 100644 index 0000000..2f93ab7 --- /dev/null +++ b/pair_correlation_calc.c @@ -0,0 +1,83 @@ +/* + * calcultae pair correlation function + * + * 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 *stat; + int slots; + int i; + double dr; + int fd; + + if(argc!=3) { + usage(argv[0]); + return -1; + } + + ret=moldyn_read_save_file(&moldyn,argv[1]); + if(ret) { + printf("[pair corr calc] exit!\n"); + return ret; + } + + dr=atof(argv[2]); + slots=(int)(moldyn.cutoff/dr); + + stat=(double *)malloc(3*slots*sizeof(double)); + if(stat==NULL) { + perror("[pair corr calc] alloc mem"); + return -1; + } + + calculate_pair_correlation(&moldyn,dr,stat); + + fd=open("pair_corr_func_ab.txt", + O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR); + dprintf(fd,"# type a - type b bonds\n"); + for(i=0;i