X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=258985eb412a83ec6d3dd466c1c76c3c3a624c2c;hb=95cfeec6fbbfa975d5ac5b99ec3f7386ca3d6071;hp=3c4d87399dcaceffb42f6d91ffe97a58bd8e841f;hpb=c40d54eb3e319b17b2f6174c4eddcfd6ee3a407b;p=physik%2Fposic.git

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;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,