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)
postproc: $(DEPS)
+pair_correlation_calc: $(DEPS)
+
.PHONY:clean
clean:
rm -vf sic postproc fluctuation_calc *.o */*.o
}
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)");
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;i<lc->cells;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;i<moldyn->count;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=(j<lc->dnlc)?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<slots;i++) {
+ /*
+ * normalization: 4 pi r r dr
+ * here: not double counting pairs -> 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;
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);
--- /dev/null
+#!/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
+
--- /dev/null
+/*
+ * calcultae pair correlation function
+ *
+ * author: frank.zirkelbach@physik.uni-augsburg.de
+ *
+ */
+
+#define _GNU_SOURCE
+#include <stdio.h>
+#include <stdlib.h>
+#include <unistd.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+
+#include "moldyn.h"
+
+int usage(char *prog) {
+
+ printf("\nusage:\n");
+ printf(" %s <save file> <dr>\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<slots;i++)
+ dprintf(fd,"%f %f\n",i*dr,stat[i]);
+ close(fd);
+
+ fd=open("pair_corr_func_aa.txt",
+ O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+ dprintf(fd,"# type a - type a bonds\n");
+ for(i=0;i<slots;i++)
+ dprintf(fd,"%f %f\n",i*dr,stat[slots+i]);
+ close(fd);
+
+ fd=open("pair_corr_func_bb.txt",
+ O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+ dprintf(fd,"# type a - type b bonds\n");
+ for(i=0;i<slots;i++)
+ dprintf(fd,"%f %f\n",i*dr,stat[2*slots+i]);
+ close(fd);
+
+ free(stat);
+
+ return 0;
+}
projects / physik / posic.git / commitdiff