X-Git-Url: https://www.hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=997c32ccb24d44607d746d1e6aed87027e9de6cd;hb=b5b47daaa3718c4dec2056fe5147668023575b8e;hp=92c4e58f3d78957a03704c40fbe62af686155a9b;hpb=fd4ecbff76b2e4fc70ca88e64fb76c85fa953888;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index 92c4e58..997c32c 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -20,6 +20,58 @@
 #include "moldyn.h"
 #include "report/report.h"
 
+/*
+ * global variables, pse and atom colors (only needed here)
+ */ 
+
+static char *pse_name[]={
+	"*",
+	"H",
+	"He",
+	"Li",
+	"Be",
+	"B",
+	"C",
+	"N",
+	"O",
+	"F",
+	"Ne",
+	"Na",
+	"Mg",
+	"Al",
+	"Si",
+	"P",
+	"S",
+	"Cl",
+	"Ar",
+};
+
+static char *pse_col[]={
+	"*",
+	"White",
+	"He",
+	"Li",
+	"Be",
+	"B",
+	"Gray",
+	"N",
+	"Blue",
+	"F",
+	"Ne",
+	"Na",
+	"Mg",
+	"Al",
+	"Yellow",
+	"P",
+	"S",
+	"Cl",
+	"Ar",
+};
+
+/*
+ * the moldyn functions
+ */
+
 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 
 	printf("[moldyn] init\n");
@@ -74,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;
@@ -312,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;
@@ -427,7 +491,6 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 		         moldyn->vlsdir);
 		system(sc);
 	}
-	if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 
 	return 0;
 }
@@ -558,6 +621,35 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
 	return 0;
 }
 
+int del_atom(t_moldyn *moldyn,int tag) {
+
+	t_atom *new,*old;
+	int cnt;
+
+	old=moldyn->atom;
+
+	new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
+	if(!new) {
+		perror("[moldyn]malloc (del atom)");
+		return -1;
+	}
+
+	for(cnt=0;cnt<tag;cnt++)
+		new[cnt]=old[cnt];
+	
+	for(cnt=tag+1;cnt<moldyn->count;cnt++) {
+		new[cnt-1]=old[cnt];
+		new[cnt-1].tag=cnt-1;
+	}
+
+	moldyn->count-=1;
+	moldyn->atom=new;
+
+	free(old);
+
+	return 0;
+}
+
 /* cubic init */
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
@@ -1371,16 +1463,13 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	moldyn->tau_square=moldyn->tau*moldyn->tau;
 	moldyn->cutoff_square=moldyn->cutoff*moldyn->cutoff;
 
-	/* energy scaling factor */
-	energy_scale=moldyn->count*EV;
-
 	/* get current time */
 	gettimeofday(&t1,NULL);
 
 	/* calculate initial forces */
 	potential_force_calc(moldyn);
 #ifdef DEBUG
-return 0;
+//return 0;
 #endif
 
 	/* some stupid checks before we actually start calculating bullshit */
@@ -1413,6 +1502,9 @@ return 0;
 		moldyn->tau_square=moldyn->tau*moldyn->tau;
 		moldyn->time_steps=sched->runs[sched->count];
 
+		/* energy scaling factor (might change!) */
+		energy_scale=moldyn->count*EV;
+
 	/* integration according to schedule */
 
 	for(i=0;i<moldyn->time_steps;i++) {
@@ -1435,7 +1527,7 @@ return 0;
 
 		/* check for log & visualization */
 		if(e) {
-			if(!(i%e))
+			if(!(moldyn->total_steps%e))
 				dprintf(moldyn->efd,
 				        "%f %f %f %f\n",
 				        moldyn->time,moldyn->ekin/energy_scale,
@@ -1443,7 +1535,7 @@ return 0;
 				        get_total_energy(moldyn)/energy_scale);
 		}
 		if(m) {
-			if(!(i%m)) {
+			if(!(moldyn->total_steps%m)) {
 				momentum=get_total_p(moldyn);
 				dprintf(moldyn->mfd,
 				        "%f %f %f %f %f\n",moldyn->time,
@@ -1452,7 +1544,7 @@ return 0;
 			}
 		}
 		if(p) {
-			if(!(i%p)) {
+			if(!(moldyn->total_steps%p)) {
 				dprintf(moldyn->pfd,
 				        "%f %f %f %f %f\n",moldyn->time,
 				         moldyn->p/BAR,moldyn->p_avg/BAR,
@@ -1460,14 +1552,14 @@ return 0;
 			}
 		}
 		if(t) {
-			if(!(i%t)) {
+			if(!(moldyn->total_steps%t)) {
 				dprintf(moldyn->tfd,
 				        "%f %f %f\n",
 				        moldyn->time,moldyn->t,moldyn->t_avg);
 			}
 		}
 		if(s) {
-			if(!(i%s)) {
+			if(!(moldyn->total_steps%s)) {
 				snprintf(dir,128,"%s/s-%07.f.save",
 				         moldyn->vlsdir,moldyn->time);
 				fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
@@ -1482,14 +1574,13 @@ return 0;
 			}	
 		}
 		if(v) {
-			if(!(i%v)) {
-				visual_atoms(&(moldyn->vis),moldyn->time,
-				             moldyn->atom,moldyn->count);
+			if(!(moldyn->total_steps%v)) {
+				visual_atoms(moldyn);
 			}
 		}
 
 		/* display progress */
-		if(!(i%10)) {
+		if(!(moldyn->total_steps%10)) {
 			/* get current time */
 			gettimeofday(&t2,NULL);
 
@@ -1820,7 +1911,13 @@ int potential_force_calc(t_moldyn *moldyn) {
 	}
 
 #ifdef DEBUG
-	printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+	//printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
+	if(moldyn->time>DSTART&&moldyn->time<DEND) {
+		printf("force:\n");
+		printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
+		printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
+		printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+	}
 #endif
 
 	/* calculate global virial */
@@ -1975,6 +2072,157 @@ int analyze_bonds(t_moldyn *moldyn) {
 	
 	
 
+	return 0;
+}
+
+/*
+ * visualization code
+ */
+
+int visual_init(t_moldyn *moldyn,char *filebase) {
+
+	strncpy(moldyn->vis.fb,filebase,128);
+
+	return 0;
+}
+
+int visual_atoms(t_moldyn *moldyn) {
+
+	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,moldyn->time);
+	fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+	if(fd<0) {
+		perror("open visual save file fd");
+		return -1;
+	}
+
+	/* write the actual data file */
+
+	// povray header
+	dprintf(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++) {
+		// 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,
+		        dim.x/2,-dim.y/2,-dim.z/2);
+		dprintf(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",
+		        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",
+		        -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",
+		        -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",
+		        -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",
+		        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",
+		        -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",
+		        -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",
+		        -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",
+		        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",
+		        dim.x/2,dim.y/2,dim.z/2,
+		        dim.x/2,dim.y/2,-dim.z/2);
+	}
+
+	close(fd);
+
 	return 0;
 }