+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[0]+=d2;
+ }
+
+ dc[2]+=d2;
+ }
+
+ dc[0]*=(1.0/(6.0*moldyn->time*a_cnt));
+ dc[1]*=(1.0/(6.0*moldyn->time*b_cnt));
+ dc[2]*=(1.0/(6.0*moldyn->time*moldyn->count));
+
+ return 0;
+}
+
+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) {
+
+ t_3dvec dist;
+ double d;
+ int s;
+ t_pcc *pcc;
+
+ /* 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 cutoff */
+ if(d>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;
+ }
+
+ return 0;
+}
+
+int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
+
+ t_pcc pcc;
+ double norm;
+ int i;
+
+ pcc.dr=dr;
+ pcc.o1=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");
+ 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) {
+ pcc.stat=(double *)ptr;
+ }
+ else {
+ pcc.stat=(double *)malloc(3*pcc.o1*sizeof(double));
+ if(pcc.stat==NULL) {
+ perror("[moldyn] pair correlation malloc");
+ return -1;
+ }
+ }
+
+ memset(pcc.stat,0,3*pcc.o1*sizeof(double));
+
+ /* process */
+ process_2b_bonds(moldyn,&pcc,calculate_pair_correlation_process);
+
+ /* normalization */
+ for(i=1;i<pcc.o1;i++) {
+ // normalization: 4 pi r^2 dr
+ // here: not double counting pairs -> 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;
+ }
+ /* */
+
+ if(ptr==NULL) {
+ /* todo: store/print pair correlation function */
+ free(pcc.stat);
+ }
+
+ return 0;
+}
+
+int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+ void *data,u8 bc) {
+
+ t_ba *ba;
+ t_3dvec dist;
+ double d;
+
+ if(itom->tag>=jtom->tag)
+ return 0;
+
+ /* 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 cutoff */
+ if(d>moldyn->cutoff_square)
+ return 0;
+
+ /* check for potential bond */
+ if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+ return 0;
+
+ /* now count this bonding ... */
+ ba=data;
+
+ /* increase total bond counter
+ */
+ ba->tcnt+=1;
+
+ 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;
+
+ return 0;
+}
+
+int bond_analyze(t_moldyn *moldyn,double *quality) {
+
+ int qcnt4;
+ int qcnt3;
+ int ccnt4;
+ int ccnt3;
+ int bcnt;
+ 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));
+
+ 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));
+
+ ba.tcnt=0;
+ qcnt4=0; qcnt3=0;
+ ccnt4=0; ccnt3=0;
+ bcnt=0;
+
+ atom=moldyn->atom;
+
+ process_2b_bonds(moldyn,&ba,bond_analyze_process);
+
+ for(i=0;i<moldyn->count;i++) {
+ if(atom[i].brand==0) {
+ if((ba.acnt[i]==0)&(ba.bcnt[i]==4))
+ qcnt4+=4;
+ if((ba.acnt[i]==0)&(ba.bcnt[i]==3))
+ qcnt3+=3;
+ }
+ else {
+ if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
+ qcnt4+=4;
+ ccnt4+=1;
+ }
+ if((ba.acnt[i]==3)&(ba.bcnt[i]==0)) {
+ qcnt3+=4;
+ ccnt3+=1;
+ }
+ bcnt+=1;
+ }
+ }
+
+ if(quality) {
+ quality[0]=1.0*ccnt4/bcnt;
+ quality[1]=1.0*ccnt3/bcnt;
+ }
+ else {
+ printf("[moldyn] bond analyze: %f %f\n",
+ 1.0*ccnt4/bcnt,1.0*ccnt3/bcnt);
+ }
+
+ return 0;
+}
+
+/*
+ * visualization code
+ */
+
+int visual_init(t_moldyn *moldyn,char *filebase) {
+
+ strncpy(moldyn->vis.fb,filebase,128);
+
+ return 0;
+}
+
+int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
+ void *data,u8 bc) {
+
+ t_vb *vb;
+
+ vb=data;
+
+ if(itom->tag>=jtom->tag)
+ return 0;
+
+ if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+ return 0;
+
+ if((itom->attr&ATOM_ATTR_VB)|(jtom->attr&ATOM_ATTR_VB))
+ dprintf(vb->fd,"# [B] %f %f %f %f %f %f\n",
+ itom->r.x,itom->r.y,itom->r.z,
+ jtom->r.x,jtom->r.y,jtom->r.z);
+
+ return 0;
+}
+
+#ifdef VISUAL_THREAD
+void *visual_atoms(void *ptr) {
+#else
+int visual_atoms(t_moldyn *moldyn) {
+#endif
+
+ int i;
+ char file[128+64];
+ t_3dvec dim;
+ double help;
+ t_visual *v;
+ t_atom *atom;
+ t_vb vb;
+ t_3dvec strain;
+#ifdef VISUAL_THREAD
+ t_moldyn *moldyn;
+
+ moldyn=ptr;
+#endif
+
+ v=&(moldyn->vis);
+ dim.x=v->dim.x;
+ dim.y=v->dim.y;
+ dim.z=v->dim.z;
+ atom=moldyn->atom;
+
+ help=(dim.x+dim.y);
+
+ sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time);
+ vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+ if(vb.fd<0) {
+ perror("open visual save file fd");
+#ifndef VISUAL_THREAD
+ return -1;
+#endif
+ }
+
+ /* write the actual data file */