+int pair_correlation_init(t_moldyn *moldyn,double dr) {
+
+
+ return 0;
+}
+
+int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc) {
+
+ int i;
+ t_atom *atom;
+ t_3dvec dist;
+ double d2;
+ int a_cnt;
+ int b_cnt;
+
+ atom=moldyn->atom;
+ dc[0]=0;
+ dc[1]=0;
+ dc[2]=0;
+ a_cnt=0;
+ b_cnt=0;
+
+ for(i=0;i<moldyn->count;i++) {
+
+ v3_sub(&dist,&(atom[i].r),&(atom[i].r_0));
+ check_per_bound(moldyn,&dist);
+ d2=v3_absolute_square(&dist);
+
+ if(atom[i].brand) {
+ b_cnt+=1;
+ dc[1]+=d2;
+ }
+ else {
+ a_cnt+=1;
+ }
+
+ dc[2]+=d2;
+ }
+
+ return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+ 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;
+
+ lc=&(moldyn->lc);
+
+ slots=moldyn->cutoff/dr;
+ o=2*slots;
+
+ if(slots*dr<=moldyn->cutoff)
+ printf("[moldyn] WARNING: pcc (low #slots)\n");
+
+ printf("[moldyn] pair correlation calc info:\n");
+ printf(" time: %f\n",moldyn->time);
+ printf(" count: %d\n",moldyn->count);
+ printf(" cutoff: %f\n",moldyn->cutoff);
+ printf(" temperature: cur=%f avg=%f\n",moldyn->t,moldyn->t_avg);
+
+ 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);
+
+ itom=moldyn->atom;
+
+ for(i=0;i<moldyn->count;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);
+
+ /* brand of atom i */
+ ibrand=itom[i].brand;
+
+ for(j=0;j<27;j++) {
+
+ bc=(j<lc->dnlc)?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
+ /* only count pairs once,
+ * skip same atoms */
+ 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 or equal cutoff */
+ if(d>=moldyn->cutoff_square)
+ continue;
+
+ /* 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;
+ }
+
+ 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
+ }
+ }
+
+ /* 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;
+}
+