+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;
+
+ d=sqrt(d);
+
+ /* now count this bonding ... */
+ ba=data;
+
+ /* increase total bond counter
+ * ... double counting!
+ */
+ ba->tcnt+=2;
+
+ 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) {
+
+ // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
+
+ int qcnt;
+ int ccnt,cset;
+ 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;
+ qcnt=0;
+ ccnt=0;
+ cset=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))
+ qcnt+=4;
+ }
+ else {
+ if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
+ qcnt+=4;
+ ccnt+=1;
+ }
+ cset+=1;
+ }
+ }
+
+ printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
+ printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
+
+ if(quality) {
+ quality[0]=1.0*ccnt/cset;
+ quality[1]=1.0*qcnt/ba.tcnt;
+ }
+ else {
+ printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
+ printf("[moldyn] bond analyze: tot_q=%f\n",1.0*qcnt/ba.tcnt);
+ }
+
+ return 0;
+}
+
+/*
+ * visualization code
+ */
+
+int visual_init(t_moldyn *moldyn,char *filebase) {
+
+ strncpy(moldyn->vis.fb,filebase,128);
+
+ return 0;
+}
+
+int visual_atoms(t_moldyn *moldyn) {
+
+ int i,j,fd;
+ char file[128+64];
+ t_3dvec dim;
+ double help;
+ t_visual *v;
+ t_atom *atom;
+ t_atom *btom;
+ t_linkcell *lc;
+#ifdef STATIC_LISTS
+ int *neighbour[27];
+ int p;
+#else
+ t_list neighbour[27];
+#endif
+ u8 bc;
+ t_3dvec dist;
+ double d2;
+ u8 brand;
+
+ v=&(moldyn->vis);
+ dim.x=v->dim.x;
+ dim.y=v->dim.y;
+ dim.z=v->dim.z;
+ atom=moldyn->atom;
+ lc=&(moldyn->lc);
+
+ help=(dim.x+dim.y);
+
+ sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
+ fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+ if(fd<0) {
+ perror("open visual save file fd");
+ return -1;
+ }
+
+ /* write the actual data file */
+
+ // povray header
+ dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+ moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
+
+ // atomic configuration
+ for(i=0;i<moldyn->count;i++) {
+ // atom type, positions, color and kinetic energy
+ dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
+ atom[i].r.x,
+ atom[i].r.y,
+ atom[i].r.z,
+ pse_col[atom[i].element],
+ atom[i].ekin);
+
+ /*
+ * bond detection should usually be done by potential
+ * functions. brrrrr! EVIL!
+ *
+ * todo: potentials need to export a 'find_bonds' function!
+ */
+
+ // bonds between atoms
+ if(!(atom[i].attr&ATOM_ATTR_VB))
+ continue;
+ link_cell_neighbour_index(moldyn,
+ (atom[i].r.x+moldyn->dim.x/2)/lc->x,
+ (atom[i].r.y+moldyn->dim.y/2)/lc->y,
+ (atom[i].r.z+moldyn->dim.z/2)/lc->z,
+ neighbour);
+ for(j=0;j<27;j++) {
+ bc=j<lc->dnlc?0:1;
+#ifdef STATIC_LISTS
+ p=0;
+ while(neighbour[j][p]!=0) {
+ btom=&(atom[neighbour[j][p]]);
+ p++;
+#else
+ list_reset_f(&neighbour[j]);
+ if(neighbour[j].start==NULL)
+ continue;
+ do {
+ btom=neighbour[j].current->data;
+#endif
+ if(btom==&atom[i]) // skip identical atoms
+ continue;
+ //if(btom<&atom[i]) // skip half of them
+ // continue;
+ v3_sub(&dist,&(atom[i].r),&(btom->r));
+ if(bc) check_per_bound(moldyn,&dist);
+ d2=v3_absolute_square(&dist);
+ brand=atom[i].brand;
+ if(brand==btom->brand) {
+ if(d2>moldyn->bondlen[brand])
+ continue;
+ }
+ else {
+ if(d2>moldyn->bondlen[2])
+ continue;
+ }
+ dprintf(fd,"# [B] %f %f %f %f %f %f\n",
+ atom[i].r.x,atom[i].r.y,atom[i].r.z,
+ btom->r.x,btom->r.y,btom->r.z);
+#ifdef STATIC_LISTS
+ }
+#else
+ } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
+#endif
+ }
+ }
+
+ // boundaries
+ if(dim.x) {
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,-dim.y/2,-dim.z/2,
+ dim.x/2,-dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,-dim.y/2,-dim.z/2,
+ -dim.x/2,dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dim.x/2,dim.y/2,-dim.z/2,
+ dim.x/2,-dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,dim.y/2,-dim.z/2,
+ dim.x/2,dim.y/2,-dim.z/2);
+
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,-dim.y/2,dim.z/2,
+ dim.x/2,-dim.y/2,dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,-dim.y/2,dim.z/2,
+ -dim.x/2,dim.y/2,dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dim.x/2,dim.y/2,dim.z/2,
+ dim.x/2,-dim.y/2,dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,dim.y/2,dim.z/2,
+ dim.x/2,dim.y/2,dim.z/2);
+
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,-dim.y/2,dim.z/2,
+ -dim.x/2,-dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ -dim.x/2,dim.y/2,dim.z/2,
+ -dim.x/2,dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dim.x/2,-dim.y/2,dim.z/2,
+ dim.x/2,-dim.y/2,-dim.z/2);
+ dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+ dim.x/2,dim.y/2,dim.z/2,
+ dim.x/2,dim.y/2,-dim.z/2);
+ }
+
+ close(fd);
+
+ return 0;
+}
+