more interstitial testing, added bond visualization

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 3c4d873..258985e 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -126,6 +126,18 @@ int set_cutoff(t_moldyn *moldyn,double cutoff) {
        return 0;
 }
 
+int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm) {
+
+       moldyn->bondlen[0]=b0*b0;
+       moldyn->bondlen[1]=b1*b1;
+       if(bm<0)
+               moldyn->bondlen[2]=b0*b1;
+       else
+               moldyn->bondlen[2]=bm*bm;
+
+       return 0;
+}
+
 int set_temperature(t_moldyn *moldyn,double t_ref) {
 
        moldyn->t_ref=t_ref;
@@ -364,7 +376,7 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
                        break;
                case VISUAL_STEP:
                        moldyn->vwrite=timer;
-                       ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
+                       ret=visual_init(moldyn,moldyn->vlsdir);
                        if(ret<0) {
                                printf("[moldyn] visual init failure\n");
                                return ret;
@@ -1563,8 +1575,7 @@ return 0;
                }
                if(v) {
                        if(!(moldyn->total_steps%v)) {
-                               visual_atoms(&(moldyn->vis),moldyn->time,
-                                            moldyn->atom,moldyn->count);
+                               visual_atoms(moldyn);
                        }
                }
 
@@ -2062,30 +2073,39 @@ int analyze_bonds(t_moldyn *moldyn) {
  * visualization code
  */
 
-int visual_init(t_visual *v,char *filebase) {
+int visual_init(t_moldyn *moldyn,char *filebase) {
 
-       char file[128+8];
-
-       strncpy(v->fb,filebase,128);
-       memset(file,0,128+8);
+       strncpy(moldyn->vis.fb,filebase,128);
 
        return 0;
 }
 
-int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
+int visual_atoms(t_moldyn *moldyn) {
 
-       int i,fd;
+       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;
+       t_list neighbour[27];
+       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,time);
+       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");
@@ -2093,15 +2113,67 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        }
 
        /* write the actual data file */
+
+       // povray header
        dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
-               n,time,help/40.0,help/40.0,-0.8*help);
-       for(i=0;i<n;i++)
+               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++) {
+                       list_reset_f(&neighbour[j]);
+                       if(neighbour[j].start==NULL)
+                               continue;
+                       bc=j<lc->dnlc?0:1;
+                       do {
+                               btom=neighbour[j].current->data;
+                               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);
+                       } while(list_next_f(&neighbour[j])!=L_NO_NEXT_ELEMENT);
+               }
+       }
+
+       // boundaries
        if(dim.x) {
                dprintf(fd,"# [D] %f %f %f %f %f %f\n",
                        -dim.x/2,-dim.y/2,-dim.z/2,