From: hackbard Date: Wed, 10 Oct 2007 09:29:35 +0000 (+0200) Subject: more interstitial testing, added bond visualization X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=95cfeec6fbbfa975d5ac5b99ec3f7386ca3d6071;p=physik%2Fposic.git more interstitial testing, added bond visualization --- diff --git a/moldyn.c b/moldyn.c index 3c4d873..258985e 100644 --- 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;icount,moldyn->time,help/40.0,help/40.0,-0.8*help); + + // atomic configuration + for(i=0;icount;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=jdnlc?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 --- 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 --- 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 --- 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 --- 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 " echo "atom size:" echo " -r " + echo " -B " echo "visualization volume:" echo " -nll (near lower left)" echo " -fur (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 {