X-Git-Url: https://hackdaworld.org/gitweb/?a=blobdiff_plain;f=moldyn.c;h=9a62b9f7c93d5142f5b15d57d099f4ba160604ce;hb=ea612b88a0588b8f46fafaebf3b37fb46c83c0cf;hp=ca6f84a65e858b8f21534cb853f7b85eb81460fc;hpb=e25ff194682ff5fac86c60701343103e74973bed;p=physik%2Fposic.git

diff --git a/moldyn.c b/moldyn.c
index ca6f84a..9a62b9f 100644
--- a/moldyn.c
+++ b/moldyn.c
@@ -246,6 +246,32 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 			dprintf(moldyn->efd,"# total momentum log file\n");
 			printf("total momentum (%d)\n",timer);
 			break;
+		case LOG_PRESSURE:
+			moldyn->pwrite=timer;
+			snprintf(filename,127,"%s/pressure",moldyn->vlsdir);
+			moldyn->pfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->pfd<0) {
+				perror("[moldyn] pressure log file\n");
+				return moldyn->pfd;
+			}
+			dprintf(moldyn->pfd,"# pressure log file\n");
+			printf("pressure (%d)\n",timer);
+			break;
+		case LOG_TEMPERATURE:
+			moldyn->twrite=timer;
+			snprintf(filename,127,"%s/temperature",moldyn->vlsdir);
+			moldyn->tfd=open(filename,
+			                 O_WRONLY|O_CREAT|O_EXCL,
+			                 S_IRUSR|S_IWUSR);
+			if(moldyn->tfd<0) {
+				perror("[moldyn] temperature log file\n");
+				return moldyn->tfd;
+			}
+			dprintf(moldyn->tfd,"# temperature log file\n");
+			printf("temperature (%d)\n",timer);
+			break;
 		case SAVE_STEP:
 			moldyn->swrite=timer;
 			printf("save file (%d)\n",timer);
@@ -268,18 +294,47 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
 				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;
+			if(moldyn->efd) {
+				snprintf(filename,127,"%s/e_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->epfd=open(filename,
+			 	                 O_WRONLY|O_CREAT|O_EXCL,
+			  	                 S_IRUSR|S_IWUSR);
+				if(moldyn->epfd<0) {
+					perror("[moldyn] energy plot fd open");
+					return moldyn->epfd;
+				}
+				dprintf(moldyn->epfd,e_plot_script);
+				close(moldyn->epfd);
+			}
+			if(moldyn->pfd) {
+				snprintf(filename,127,"%s/pressure_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->ppfd=open(filename,
+			 	                  O_WRONLY|O_CREAT|O_EXCL,
+			  	                  S_IRUSR|S_IWUSR);
+				if(moldyn->ppfd<0) {
+					perror("[moldyn] p plot fd open");
+					return moldyn->ppfd;
+				}
+				dprintf(moldyn->ppfd,pressure_plot_script);
+				close(moldyn->ppfd);
+			}
+			if(moldyn->tfd) {
+				snprintf(filename,127,"%s/temperature_plot.scr",
+				         moldyn->vlsdir);
+				moldyn->tpfd=open(filename,
+			 	                  O_WRONLY|O_CREAT|O_EXCL,
+			  	                  S_IRUSR|S_IWUSR);
+				if(moldyn->tpfd<0) {
+					perror("[moldyn] t plot fd open");
+					return moldyn->tpfd;
+				}
+				dprintf(moldyn->tpfd,temperature_plot_script);
+				close(moldyn->tpfd);
 			}
 			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);
@@ -294,13 +349,35 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
 	char sc[256];
 
 	printf("[moldyn] log shutdown\n");
-	if(moldyn->efd) close(moldyn->efd);
+	if(moldyn->efd) {
+		close(moldyn->efd);
+		if(moldyn->rfd) {
+			dprintf(moldyn->rfd,report_energy);
+			snprintf(sc,255,"cd %s && gnuplot e_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+		}
+	}
 	if(moldyn->mfd) close(moldyn->mfd);
+	if(moldyn->pfd) {
+		close(moldyn->pfd);
+		if(moldyn->rfd)
+			dprintf(moldyn->rfd,report_pressure);
+			snprintf(sc,255,"cd %s && gnuplot pressure_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+	}
+	if(moldyn->tfd) {
+		close(moldyn->tfd);
+		if(moldyn->rfd)
+			dprintf(moldyn->rfd,report_temperature);
+			snprintf(sc,255,"cd %s && gnuplot temperature_plot.scr",
+			         moldyn->vlsdir);
+			system(sc);
+	}
 	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);
@@ -676,9 +753,11 @@ double pressure_calc(t_moldyn *moldyn) {
 	t_virial *virial;
 
 	/*
-	 * P = 1/(3V) sum_i ( p_i^2 / 2m + f_i r_i )
-	 *
-	 * virial = f_i r_i
+	 * PV = NkT + <W>
+	 * W = 1/3 sum_i f_i r_i
+	 * virial = sum_i f_i r_i
+	 * 
+	 * => P = (2 Ekin + virial) / (3V)
 	 */
 
 	v=0.0;
@@ -698,9 +777,20 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 
 	t_3dvec dim,*tp;
 	double u,p;
-	double scale;
+	double scale,dv;
 	t_atom *store;
 
+	/*
+	 * dU = - p dV
+	 *
+	 * => p = - dU/dV
+	 *
+	 * dV: dx,y,z = 0.001 x,y,z
+	 */
+
+	scale=1.00001;
+printf("\n\nP-DEBUG:\n");
+
 	tp=&(moldyn->tp);
 	store=malloc(moldyn->count*sizeof(t_atom));
 	if(store==NULL) {
@@ -714,27 +804,28 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	dim=moldyn->dim;
 
 	/* derivative with respect to x direction */
-	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);
+	dv=0.00001*moldyn->dim.x*moldyn->dim.y*moldyn->dim.z;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->x=(moldyn->energy-u)/moldyn->dv;
+	tp->x=(moldyn->energy-u)/dv;
 	p=tp->x*tp->x;
+printf("e: %f eV de: %f eV dV: %f A^3\n",moldyn->energy/moldyn->count/EV,(moldyn->energy-u)/moldyn->count/EV,dv);
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
 	/* derivative with respect to y direction */
-	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);
+	dv=0.00001*moldyn->dim.y*moldyn->dim.x*moldyn->dim.z;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->y=(moldyn->energy-u)/moldyn->dv;
+	tp->y=(moldyn->energy-u)/dv;
 	p+=tp->y*tp->y;
 
 	/* restore atomic configuration + dim */
@@ -742,37 +833,19 @@ double thermodynamic_pressure_calc(t_moldyn *moldyn) {
 	moldyn->dim=dim;
 
 	/* derivative with respect to z direction */
-	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);
+	dv=0.00001*moldyn->dim.z*moldyn->dim.x*moldyn->dim.y;
 	link_cell_shutdown(moldyn);
 	link_cell_init(moldyn,QUIET);
 	potential_force_calc(moldyn);
-	tp->z=(moldyn->energy-u)/moldyn->dv;
+	tp->z=(moldyn->energy-u)/dv;
 	p+=tp->z*tp->z;
 
 	/* restore atomic configuration + dim */
 	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
 	moldyn->dim=dim;
 
-	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_atoms(moldyn,scale,1,1,1);
-	link_cell_shutdown(moldyn);
-	link_cell_init(moldyn,QUIET);
-	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);
-
-	/* restore atomic configuration + dim */
-	memcpy(moldyn->atom,store,moldyn->count*sizeof(t_atom));
-	moldyn->dim=dim;
-
 	/* restore energy */
 	moldyn->energy=u;
 
@@ -1109,14 +1182,15 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
 int moldyn_integrate(t_moldyn *moldyn) {
 
 	int i;
-	unsigned int e,m,s,v;
-	t_3dvec p;
+	unsigned int e,m,s,v,p,t;
+	t_3dvec momentum;
 	t_moldyn_schedule *sched;
 	t_atom *atom;
 	int fd;
 	char dir[128];
 	double ds;
 	double energy_scale;
+	//double tp;
 
 	sched=&(moldyn->schedule);
 	atom=moldyn->atom;
@@ -1129,6 +1203,8 @@ int moldyn_integrate(t_moldyn *moldyn) {
 	m=moldyn->mwrite;
 	s=moldyn->swrite;
 	v=moldyn->vwrite;
+	p=moldyn->pwrite;
+	t=moldyn->twrite;
 
 	/* sqaure of some variables */
 	moldyn->tau_square=moldyn->tau*moldyn->tau;
@@ -1179,7 +1255,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		update_e_kin(moldyn);
 		temperature_calc(moldyn);
 		pressure_calc(moldyn);
-		//thermodynamic_pressure_calc(moldyn);
+		//tp=thermodynamic_pressure_calc(moldyn);
 
 		/* p/t scaling */
 		if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
@@ -1198,9 +1274,23 @@ int moldyn_integrate(t_moldyn *moldyn) {
 		}
 		if(m) {
 			if(!(i%m)) {
-				p=get_total_p(moldyn);
+				momentum=get_total_p(moldyn);
 				dprintf(moldyn->mfd,
-				        "%f %f\n",moldyn->time,v3_norm(&p));
+				        "%f %f %f %f %f\n",moldyn->time,
+				        momentum.x,momentum.y,momentum.z,
+				        v3_norm(&momentum));
+			}
+		}
+		if(p) {
+			if(!(i%p)) {
+				dprintf(moldyn->pfd,
+				        "%f %f\n",moldyn->time,moldyn->p/ATM);
+			}
+		}
+		if(t) {
+			if(!(i%t)) {
+				dprintf(moldyn->tfd,
+				        "%f %f\n",moldyn->time,moldyn->t);
 			}
 		}
 		if(s) {
@@ -1221,13 +1311,17 @@ int moldyn_integrate(t_moldyn *moldyn) {
 			if(!(i%v)) {
 				visual_atoms(&(moldyn->vis),moldyn->time,
 				             moldyn->atom,moldyn->count);
-				printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f",
-				       sched->count,i,
-				       moldyn->t,moldyn->p/ATM,moldyn->volume);
-				fflush(stdout);
 			}
 		}
 
+		/* display progress */
+		if(!(i%10)) {
+			printf("\rsched: %d, steps: %d, T: %f, P: %f V: %f",
+			       sched->count,i,
+			       moldyn->t,moldyn->p/ATM,moldyn->volume);
+			fflush(stdout);
+		}
+
 		/* increase absolute time */
 		moldyn->time+=moldyn->tau;