mdrun iinit. untested + not working by now!
[physik/posic.git] / moldyn.c
index a1bfcae..b52d51d 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -131,24 +131,13 @@ int set_int_alg(t_moldyn *moldyn,u8 algo) {
 int set_cutoff(t_moldyn *moldyn,double cutoff) {
 
        moldyn->cutoff=cutoff;
+       moldyn->cutoff_square=cutoff*cutoff;
 
        printf("[moldyn] cutoff [A]: %f\n",moldyn->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;
@@ -1581,7 +1570,6 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
        /* sqaure of some variables */
        moldyn->tau_square=moldyn->tau*moldyn->tau;
-       moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 
        /* get current time */
        gettimeofday(&t1,NULL);
@@ -2414,7 +2402,8 @@ int get_line(int fd,char *line,int max) {
                ret=read(fd,line+count,1);
                if(ret<=0) return ret;
                if(line[count]=='\n') {
-                       line[count]='\0';
+                       memset(line+count,0,max-count-1);
+                       //line[count]='\0';
                        return count+1;
                }
                count+=1;
@@ -2705,39 +2694,48 @@ int visual_init(t_moldyn *moldyn,char *filebase) {
        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;
+}
+
 int visual_atoms(t_moldyn *moldyn) {
 
-       int i,j,fd;
+       int i;
        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;
+       t_vb vb;
 
        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) {
+       vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       if(vb.fd<0) {
                perror("open visual save file fd");
                return -1;
        }
@@ -2745,118 +2743,65 @@ int visual_atoms(t_moldyn *moldyn) {
        /* write the actual data file */
 
        // povray header
-       dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+       dprintf(vb.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++) {
+       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
-               }
-       }
-
+               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);
+       
+       // bonds between atoms
+       process_2b_bonds(moldyn,&vb,visual_bonds_process);
+       
        // boundaries
        if(dim.x) {
-               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd,"# [D] %f %f %f %f %f %f\n",
+               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(fd);
+       close(vb.fd);
 
        return 0;
 }