X-Git-Url: https://www.hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.c;h=817a72788b45f2572258873aaa8ca64cace0397f;hp=94d596e91878237cf41a127cd8aa947cbcd18997;hb=f3e4193447ac49a8953515910d4b4e6ce2c7608b;hpb=fb951c04e522e4637618bf622fc67194c2a7b15f

diff --git a/moldyn.c b/moldyn.c
index 94d596e..817a727 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -121,7 +121,7 @@ int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize) {
 		moldyn->vis.dim.z=z;
 	}
 
-	moldyn->dv=0.0001*moldyn->volume;
+	moldyn->dv=0.000001*moldyn->volume;
 
 	printf("[moldyn] dimensions in A and A^3 respectively:\n");
 	printf("  x: %f\n",moldyn->dim.x);
@@ -199,6 +199,14 @@ int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
 
 	return 0;
 }
+
+int moldyn_set_report(t_moldyn *moldyn,char *author,char *title) {
+
+	strncpy(moldyn->rauthor,author,63);
+	strncpy(moldyn->rtitle,title,63);
+
+	return 0;
+}
 	
 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 
@@ -247,6 +255,28 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 			}
 			printf("visual file (%d)\n",timer);
 			break;
+		case CREATE_REPORT:
+			snprintf(filename,127,"%s/report.tex",moldyn->vlsdir);
+			moldyn->rfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->rfd<0) {
+				perror("[moldyn] report fd open");	
+				return moldyn->rfd;
+			}
+			snprintf(filename,127,"%s/plot.scr",moldyn->vlsdir);
+			moldyn->pfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->pfd<0) {
+				perror("[moldyn] plot fd open");
+				return moldyn->pfd;
+			}
+			dprintf(moldyn->rfd,report_start,
+			        moldyn->rauthor,moldyn->rtitle);
+			dprintf(moldyn->pfd,plot_script);
+			close(moldyn->pfd);
+			break;
 		default:
 			printf("unknown log type: %02x\n",type);
 			return -1;
@@ -257,9 +287,23 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 
 int moldyn_log_shutdown(t_moldyn *moldyn) {
 
+	char sc[256];
+
 	printf("[moldyn] log shutdown\n");
 	if(moldyn->efd) close(moldyn->efd);
 	if(moldyn->mfd) close(moldyn->mfd);
+	if(moldyn->rfd) {
+		dprintf(moldyn->rfd,report_end);
+		close(moldyn->rfd);
+		snprintf(sc,255,"cd %s && gnuplot plot.scr",moldyn->vlsdir);
+		system(sc);
+		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+		system(sc);
+		snprintf(sc,255,"cd %s && pdflatex report",moldyn->vlsdir);
+		system(sc);
+		snprintf(sc,255,"cd %s && dvipdf report",moldyn->vlsdir);
+		system(sc);
+	}
 	if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 
 	return 0;
@@ -299,10 +343,10 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
 	switch(type) {
 		case CUBIC:
-			ret=cubic_init(a,b,c,atom,&origin);
+			ret=cubic_init(a,b,c,lc,atom,NULL);
 			break;
 		case FCC:
-			ret=fcc_init(a,b,c,lc,atom,&origin);
+			ret=fcc_init(a,b,c,lc,atom,NULL);
 			break;
 		case DIAMOND:
 			ret=diamond_init(a,b,c,lc,atom,&origin);
@@ -340,9 +384,38 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
 	int count;
-	t_3dvec r
+	t_3dvec r;
+	int i,j,k;
+	t_3dvec o;
+
+	count=0;
+	if(origin)
+		v3_copy(&o,origin);
+	else
+		v3_zero(&o);
 
-HIER WEITER !
+	r.x=o.x;
+	for(i=0;i<a;i++) {
+		r.y=o.y;
+		for(j=0;j<b;j++) {
+			for(k=0;k<c;k++) {
+				r.z=o.z;
+				v3_copy(&(atom[count].r),&r);
+				count+=1;
+				r.z+=lc;
+			}
+			r.y+=lc;
+		}
+		r.x+=lc;
+	}
+
+	for(i=0;i<count;i++) {
+		atom[i].r.x-=(a*lc)/2.0;
+		atom[i].r.y-=(b*lc)/2.0;
+		atom[i].r.z-=(c*lc)/2.0;
+	}
+
+	return count;
 }
 
 /* fcc lattice init */
@@ -599,7 +672,6 @@ double ideal_gas_law_pressure(t_moldyn *moldyn) {
 	double p;
 
 	p=moldyn->count*moldyn->t*K_BOLTZMANN/moldyn->volume;
-printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM);
 
 	return p;
 }
@@ -607,8 +679,7 @@ printf("temp = %f => ideal gas law pressure: %f\n",moldyn->t,p/ATM);
 double pressure_calc(t_moldyn *moldyn) {
 
 	int i;
-	double p1,p;
-	double v,h;
+	double v;
 	t_virial *virial;
 
 	v=0.0;
@@ -616,18 +687,11 @@ double pressure_calc(t_moldyn *moldyn) {
 		virial=&(moldyn->atom[i].virial);
 		v+=(virial->xx+virial->yy+virial->zz);
 	}
+	v*=ONE_THIRD;
+printf("kieck mal: %f %f %f\n",v,moldyn->count*K_BOLTZMANN*moldyn->t,v/moldyn->count);
 
-	h=moldyn->count*K_BOLTZMANN*moldyn->t;
-
-	p=(h-ONE_THIRD*v);
-	p/=moldyn->volume;
-
-	p1=(moldyn->count*K_BOLTZMANN*moldyn->t-ONE_THIRD*moldyn->vt1);
-	p1/=moldyn->volume;
-
-	printf("debug: vt1=%f v=%f nkt=%f\n",moldyn->vt1,v,h);
-
-	printf("compare pressures: %f %f %f\n",p1/ATM,p/ATM,h/moldyn->volume/ATM);
+	moldyn->p=moldyn->count*K_BOLTZMANN*moldyn->t+v;
+	moldyn->p/=moldyn->volume;
 
 	return moldyn->p;
 }	
@@ -655,6 +719,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	scale=1.0+moldyn->dv/(moldyn->dim.y*moldyn->dim.z);
 	scale_dim(moldyn,scale,TRUE,0,0);
 	scale_atoms(moldyn,scale,TRUE,0,0);
+	link_cell_shutdown(moldyn);
+	link_cell_init(moldyn);
 	potential_force_calc(moldyn);
 	tp->x=(moldyn->energy-u)/moldyn->dv;
 	p=tp->x*tp->x;
@@ -667,6 +733,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.z);
 	scale_dim(moldyn,scale,0,TRUE,0);
 	scale_atoms(moldyn,scale,0,TRUE,0);
+	link_cell_shutdown(moldyn);
+	link_cell_init(moldyn);
 	potential_force_calc(moldyn);
 	tp->y=(moldyn->energy-u)/moldyn->dv;
 	p+=tp->y*tp->y;
@@ -679,6 +747,8 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	scale=1.0+moldyn->dv/(moldyn->dim.x*moldyn->dim.y);
 	scale_dim(moldyn,scale,0,0,TRUE);
 	scale_atoms(moldyn,scale,0,0,TRUE);
+	link_cell_shutdown(moldyn);
+	link_cell_init(moldyn);
 	potential_force_calc(moldyn);
 	tp->z=(moldyn->energy-u)/moldyn->dv;
 	p+=tp->z*tp->z;
@@ -687,15 +757,19 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
-	printf("dU/dV komp addiert = %f\n",(tp->x+tp->y+tp->z)/ATM);
+	printf("dU/dV komp addiert = %f %f %f\n",tp->x,tp->y,tp->z);
 
 	scale=1.0+pow(moldyn->dv/moldyn->volume,ONE_THIRD);
 
+printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
 	scale_dim(moldyn,scale,1,1,1);
-	scale_dim(moldyn,scale,1,1,1);
+	scale_atoms(moldyn,scale,1,1,1);
+	link_cell_shutdown(moldyn);
+	link_cell_init(moldyn);
 	potential_force_calc(moldyn);
+printf("debug: %f %f\n",moldyn->atom[0].r.x,moldyn->dim.x);
 
-	printf("dU/dV einfach = %f\n",(moldyn->energy-u)/moldyn->dv/ATM);
+	printf("dU/dV einfach = %f\n",((moldyn->energy-u)/moldyn->dv)/ATM);
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
@@ -704,6 +778,9 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	/* restore energy */
 	moldyn->energy=u;
 
+	link_cell_shutdown(moldyn);
+	link_cell_init(moldyn);
+
 	return sqrt(p);
 }
 
@@ -1054,6 +1131,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	int fd;
 	char dir[128];
 	double ds;
+	double energy_scale;
 
 	sched=&(moldyn->schedule);
 	atom=moldyn->atom;
@@ -1071,6 +1149,9 @@ 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;
+
 	/* calculate initial forces */
 	potential_force_calc(moldyn);
 
@@ -1115,7 +1196,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		if(moldyn->pt_scale&(P_SCALE_BERENDSEN|P_SCALE_DIRECT))
 			scale_volume(moldyn);
 
-printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM);
+		temperature_calc(moldyn);
+		pressure_calc(moldyn);
+		//thermodynamic_pressure_calc(moldyn);
 
 		/* check for log & visualization */
 		if(e) {
@@ -1123,9 +1206,9 @@ printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM);
 				update_e_kin(moldyn);
 				dprintf(moldyn->efd,
 				        "%f %f %f %f\n",
-				        moldyn->time,moldyn->ekin,
-				        moldyn->energy,
-				        get_total_energy(moldyn));
+				        moldyn->time,moldyn->ekin/energy_scale,
+				        moldyn->energy/energy_scale,
+				        get_total_energy(moldyn)/energy_scale);
 		}
 		if(m) {
 			if(!(i%m)) {
@@ -1152,8 +1235,8 @@ printf("-> %f\n",thermodynamic_pressure_calc(moldyn)/ATM);
 			if(!(i%v)) {
 				visual_atoms(&(moldyn->vis),moldyn->time,
 				             moldyn->atom,moldyn->count);
-				printf("\rsched: %d, steps: %d, debug: %d",
-				       sched->count,i,moldyn->debug);
+				printf("\rsched: %d, steps: %d, debug: %f | %f",
+				       sched->count,i,moldyn->p/ATM,moldyn->debug/ATM);
 				fflush(stdout);
 			}
 		}
@@ -1246,9 +1329,6 @@ int potential_force_calc(t_moldyn *moldyn) {
 	/* reset energy */
 	moldyn->energy=0.0;
 
-	/* reset virial */
-	moldyn->vt1=0.0;
-
 	/* reset force, site energy and virial of every atom */
 	for(i=0;i<count;i++) {
 
@@ -1370,6 +1450,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 		}
 
 	}
+
 #ifdef DEBUG
 printf("\n\n");
 #endif
@@ -1377,10 +1458,6 @@ printf("\n\n");
 printf("\n\n");
 #endif
 
-temperature_calc(moldyn);
-pressure_calc(moldyn);
-ideal_gas_law_pressure(moldyn);
-
 	return 0;
 }
 
@@ -1390,12 +1467,12 @@ ideal_gas_law_pressure(moldyn);
 
 inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) {
 
-	a->virial.xx-=f->x*d->x;
-	a->virial.yy-=f->y*d->y;
-	a->virial.zz-=f->z*d->z;
-	a->virial.xy-=f->x*d->y;
-	a->virial.xz-=f->x*d->z;
-	a->virial.yz-=f->y*d->z;
+	a->virial.xx+=f->x*d->x;
+	a->virial.yy+=f->y*d->y;
+	a->virial.zz+=f->z*d->z;
+	a->virial.xy+=f->x*d->y;
+	a->virial.xz+=f->x*d->z;
+	a->virial.yz+=f->y*d->z;
 
 	return 0;
 }
@@ -1507,7 +1584,6 @@ int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 		v3_add(&(ai->f),&(ai->f),&force);
 		virial_calc(ai,&force,&distance);
 		virial_calc(aj,&force,&distance); /* f and d signe switched */
-		moldyn->vt1-=v3_scalar_product(&force,&distance);
 	}
 
 	return 0;