more interstitial testing, added bond visualization

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,

diff --git a/moldyn.h b/moldyn.h
index 83414dd..2d4ca48 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -49,6 +49,8 @@ typedef struct s_atom {
 
 #define ATOM_ATTR_FP   0x01    /* fixed position (bulk material) */
 #define ATOM_ATTR_HB   0x02    /* coupled to heat bath (velocity scaling) */
+#define ATOM_ATTR_VA   0x04    /* visualize this atom */
+#define ATOM_ATTR_VB   0x08    /* visualize the bond of this atom */
 
 #define ATOM_ATTR_1BP          0x10    /* single paricle potential */
 #define ATOM_ATTR_2BP          0x20    /* pair potential */
@@ -109,6 +111,7 @@ typedef struct s_moldyn {
        double cutoff;          /* cutoff radius */
        double cutoff_square;   /* square of the cutoff radius */
        double nnd;             /* nearest neighbour distance (optional) */
+       double bondlen[3];      /* bond lengthes (only 2 atomic systems) */
 
        t_linkcell lc;          /* linked cell list interface */
 
@@ -395,6 +398,7 @@ int moldyn_shutdown(t_moldyn *moldyn);
 
 int set_int_alg(t_moldyn *moldyn,u8 algo);
 int set_cutoff(t_moldyn *moldyn,double cutoff);
+int set_bondlen(t_moldyn *moldyn,double b0,double b1,double bm);
 int set_temperature(t_moldyn *moldyn,double t_ref);
 int set_pressure(t_moldyn *moldyn,double p_ref);
 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
@@ -473,8 +477,8 @@ int moldyn_bc_check(t_moldyn *moldyn);
 
 int get_line(int fd,char *line,int max);
 
-int visual_init(t_visual *v,char *filebase);
-int visual_atoms(t_visual *v,double time,t_atom *atom,int n);
+int visual_init(t_moldyn *moldyn,char *filebase);
+int visual_atoms(t_moldyn *moldyn);
 
 #endif
 

diff --git a/run b/run
index c7d454c..bd72bc2 100755 (executable)
--- a/run
+++ b/run
@@ -4,7 +4,8 @@
 if [ "$?" == "0" ]; then
        #./perms
        if [ "$1" ] ; then
-               ./visualize -d $1
+               #./visualize -w 640 -h 480 -d $1
+               ./visualize -w 640 -h 480 -d $1 -nll -2.4 -2.4 -2.4 -fur 3.8 3.8 3.8 -b -2.03 -2.03 -2.03 3.39 3.39 3.39 -r 0.6 -c 1.5 -15.0 1.5
                #rasmol -32 -nodisplay < $1/visualize.scr > /dev/null 2>&1
                ./ppm2avi $1
        fi

diff --git a/sic.c b/sic.c
index 6e6369f..20c80d1 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -138,7 +138,7 @@ int hook_add_atom(void *moldyn,void *hook_params) {
 #else
                add_atom(md,SI,M_SI,0,
 #endif
-                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+                        ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB|ATOM_ATTR_VB,
                         &r,&v);
        }
        hp->a_count+=NR_ATOMS;
@@ -198,12 +198,14 @@ int main(int argc,char **argv) {
        set_potential_params(&md,&tp);
 #endif
 
-       /* cutoff radius */
+       /* cutoff radius & bondlen */
 #ifdef ALBE
        set_cutoff(&md,ALBE_S_SI);
+       set_bondlen(&md,ALBE_S_SI,ALBE_S_C,ALBE_S_SIC);
        //set_cutoff(&md,ALBE_S_C);
 #else
        set_cutoff(&md,TM_S_SI);
+       set_bondlen(&md,TM_S_SI,TM_S_C,-1.0);
        //set_cutoff(&md,TM_S_C);
 #endif
 

diff --git a/visualize b/visualize
index b5d8f74..ab6e374 100755 (executable)
--- a/visualize
+++ b/visualize
@@ -14,6 +14,14 @@ pigment { color White }
 }
 EOF
 }
+draw_bond() {
+       cat >> temp.pov <<-EOF
+cylinder {
+<$1, $3, $2>, <$4, $6, $5>, $7
+pigment { color Blue }
+}
+EOF
+}
 
 directory="doesnt_exist____for_sure"
 width="640"
@@ -23,8 +31,10 @@ x0=""; y0=""; z0="";
 x1=""; y1=""; z1="";
 cx=""; cy=""; cz="";
 lx="0"; ly="-100"; lz="100";
+ortographic=""
 bx0=""; by0=""; bz0="";
 bx1=""; by1=""; bz1="";
+bcr="";
 
 while [ "$1" ]; do
        case "$1" in
@@ -39,6 +49,7 @@ while [ "$1" ]; do
                -o)             ortographic=1;          shift 1;;
                -b)             bx0=$2; by0=$3; bz0=$4;
                                bx1=$5; by1=$6; bz1=$7; shift 7;;
+               -B)             bcr=$2;                 shift 2;;
                *)
                                echo "options:"
                                echo "########"
@@ -49,6 +60,7 @@ while [ "$1" ]; do
                                echo "  -h <height>"
                                echo "atom size:"
                                echo "  -r <radius>"
+                               echo "  -B <bond cylinder radius>"
                                echo "visualization volume:"
                                echo "  -nll <x> <y> <z> (near lower left)"
                                echo "  -fur <x> <y> <z> (far upper right)"
@@ -125,14 +137,11 @@ EOF
        if [ -z "$bx0" ]; then
 
        if [ -z "$x0" ]; then
+
        cat $file | grep '# \[D\]' | while read foo bar x1 y1 z1 x2 y2 z2 ; do
-               cat >> temp.pov <<-EOF
-cylinder {
-<$x1, $z1, $y1>, <$x2, $z2, $y2>, 0.05
-pigment { color White }
-}
-EOF
+               draw_cyl $x1 $z1 $y1 $x2 $z2 $y2 0.05
        done
+
        else
                # manually drawing the 3x4 boundaries ...
                draw_cyl $x0 $y0 $z0 $x1 $y0 $z0
@@ -171,6 +180,38 @@ EOF
 
        fi      
 
+       # bonds
+       if [ -n "$bcr" ]; then
+
+               if [ -z "$x0" ]; then
+
+       cat $file | grep '# \[B\]' | while read foo bar x1 y1 z1 x2 y2 z2 ; do
+               draw_bond $x1 $z1 $y1 $x2 $z2 $y2 $bcr
+       done
+
+               else
+
+               export x0 y0 z0 x1 y1 z1 bcr
+               cat $file | grep '# \[B\]' | awk '\
+               BEGIN {
+                       x0=ENVIRON["x0"]; y0=ENVIRON["y0"]; z0=ENVIRON["z0"];
+                       x1=ENVIRON["x1"]; y1=ENVIRON["y1"]; z1=ENVIRON["z1"];
+                       bcr=ENVIRON["bcr"];
+               }
+               {
+                       if(($3>=x0)&&($4>=y0)&&($5>=z0)&&\
+                          ($3<=x1)&&($4<=y1)&&($5<=z1)) {
+                               print "cylinder {";
+                               print "<"$3","$5","$4">,";
+                               print "<"$6","$8","$7">, "bcr;
+                               print "pigment { color Blue }";
+                               print "}";
+                       }
+               }' >> temp.pov
+
+               fi
+       fi      
+
        # add camera and light source
        cat >> temp.pov <<-EOF
 camera {