+
+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 */
+
+ // povray header
+ dprintf(vb.fd,"# [P] %d %08.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++) {
+ v3_sub(&strain,&(atom[i].r),&(atom[i].r_0));
+ check_per_bound(moldyn,&strain);
+ // atom type, positions, color and kinetic energy
+ dprintf(vb.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);
+ sqrt(v3_absolute_square(&strain)));
+ }
+
+ // bonds between atoms
+#ifndef VISUAL_THREAD
+ process_2b_bonds(moldyn,&vb,visual_bonds_process);
+#endif
+
+ // boundaries
+ if(dim.x) {
+ dprintf(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.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(vb.fd);
+
+#ifdef VISUAL_THREAD
+ pthread_exit(NULL);
+
+}
+#else
+
+ return 0;
+}
+#endif
+
+/*
+ * fpu cntrol functions
+ */
+
+// set rounding to double (eliminates -ffloat-store!)
+int fpu_set_rtd(void) {
+
+ fpu_control_t ctrl;
+
+ _FPU_GETCW(ctrl);
+
+ ctrl&=~_FPU_EXTENDED;
+ ctrl|=_FPU_DOUBLE;
+
+ _FPU_SETCW(ctrl);
+
+ return 0;
+}
+