From: hackbard <hackbard>
Date: Wed, 14 Mar 2007 15:12:05 +0000 (+0000)
Subject: improved log/report subsystem, playing around w/ pressure, sic hook
X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=ea612b88a0588b8f46fafaebf3b37fb46c83c0cf;p=physik%2Fposic.git

improved log/report subsystem, playing around w/ pressure, sic hook
testing
---

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;
 
diff --git a/moldyn.h b/moldyn.h
index b140372..85a4950 100644
--- a/moldyn.h
+++ b/moldyn.h
@@ -133,12 +133,18 @@ typedef struct s_moldyn {
 	int efd;		/* fd for energy log */
 	unsigned int mwrite;	/* how often to log momentum */
 	int mfd;		/* fd for momentum log */
+	unsigned int pwrite;	/* how often to log pressure */
+	int pfd;		/* fd for pressure log */
+	unsigned int twrite;	/* how often to log temperature */
+	int tfd;		/* fd for temperature log */
 	unsigned int vwrite;	/* how often to visualize atom information */
 	unsigned int swrite;	/* how often to create a save file */
 	int rfd;		/* report file descriptor */
 	char rtitle[64];	/* report title */
 	char rauthor[64];	/* report author */
-	int pfd;		/* gnuplot script file descriptor */
+	int epfd;		/* energy gnuplot script file descriptor */
+	int ppfd;		/* pressure gnuplot script file descriptor */
+	int tpfd;		/* temperature gnuplot script file descriptor */
 
 	u8 status;		/* general moldyn properties */
 
@@ -201,9 +207,11 @@ typedef struct s_moldyn {
 
 #define LOG_TOTAL_ENERGY		0x01
 #define LOG_TOTAL_MOMENTUM		0x02
-#define SAVE_STEP			0x04
-#define VISUAL_STEP			0x08
-#define CREATE_REPORT			0x10
+#define LOG_PRESSURE			0x04
+#define LOG_TEMPERATURE			0x08
+#define SAVE_STEP			0x10
+#define VISUAL_STEP			0x20
+#define CREATE_REPORT			0x40
 
 #define TRUE				1
 #define FALSE				0
diff --git a/potentials/tersoff.c b/potentials/tersoff.c
index 01364ea..16f771c 100644
--- a/potentials/tersoff.c
+++ b/potentials/tersoff.c
@@ -211,12 +211,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* virial */
-	ai->virial.xx-=force.x*dist_ij.x;
-	ai->virial.yy-=force.y*dist_ij.y;
-	ai->virial.zz-=force.z*dist_ij.z;
-	ai->virial.xy-=force.x*dist_ij.y;
-	ai->virial.xz-=force.x*dist_ij.z;
-	ai->virial.yz-=force.y*dist_ij.z;
+	virial_calc(ai,&force,&dist_ij);
+	//ai->virial.xx-=force.x*dist_ij.x;
+	//ai->virial.yy-=force.y*dist_ij.y;
+	//ai->virial.zz-=force.z*dist_ij.z;
+	//ai->virial.xy-=force.x*dist_ij.y;
+	//ai->virial.xz-=force.x*dist_ij.z;
+	//ai->virial.yz-=force.y*dist_ij.z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
@@ -326,12 +327,13 @@ int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
 	v3_add(&(ai->f),&(ai->f),&force);
 
 	/* virial */
-	ai->virial.xx-=force.x*dist_ij->x;
-	ai->virial.yy-=force.y*dist_ij->y;
-	ai->virial.zz-=force.z*dist_ij->z;
-	ai->virial.xy-=force.x*dist_ij->y;
-	ai->virial.xz-=force.x*dist_ij->z;
-	ai->virial.yz-=force.y*dist_ij->z;
+	virial_calc(ai,&force,dist_ij);
+	//ai->virial.xx-=force.x*dist_ij->x;
+	//ai->virial.yy-=force.y*dist_ij->y;
+	//ai->virial.zz-=force.z*dist_ij->z;
+	//ai->virial.xy-=force.x*dist_ij->y;
+	//ai->virial.xz-=force.x*dist_ij->z;
+	//ai->virial.yz-=force.y*dist_ij->z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
@@ -380,12 +382,13 @@ if(ai==&(moldyn->atom[0])) {
 
 	/* virial - plus sign, as dist_ij = - dist_ji - (really??) */
 // TEST ... with a minus instead
-	ai->virial.xx-=force.x*dist_ij->x;
-	ai->virial.yy-=force.y*dist_ij->y;
-	ai->virial.zz-=force.z*dist_ij->z;
-	ai->virial.xy-=force.x*dist_ij->y;
-	ai->virial.xz-=force.x*dist_ij->z;
-	ai->virial.yz-=force.y*dist_ij->z;
+	virial_calc(ai,&force,dist_ij);
+	//ai->virial.xx-=force.x*dist_ij->x;
+	//ai->virial.yy-=force.y*dist_ij->y;
+	//ai->virial.zz-=force.z*dist_ij->z;
+	//ai->virial.xy-=force.x*dist_ij->y;
+	//ai->virial.xz-=force.x*dist_ij->z;
+	//ai->virial.yz-=force.y*dist_ij->z;
 
 #ifdef DEBUG
 if(ai==&(moldyn->atom[0])) {
diff --git a/report/report.h b/report/report.h
index e46ff2a..4fe42f7 100644
--- a/report/report.h
+++ b/report/report.h
@@ -34,17 +34,38 @@ static char report_start[]="\
 \\maketitle\n\
 ";
 
-static char report_end[]="\
+static char report_energy[]="\
 \\begin{figure}[!h]\n\
 \\begin{center}\n\
 \\includegraphics[width=16cm]{energy.eps}\n\
 \\caption{Kinetic, potential and total energy over time}\n\
 \\end{center}\n\
 \\end{figure}\n\
+";
+
+static char report_pressure[]="\
+\\begin{figure}[!h]\n\
+\\begin{center}\n\
+\\includegraphics[width=16cm]{pressure.eps}\n\
+\\caption{Pressure over time}\n\
+\\end{center}\n\
+\\end{figure}\n\
+";
+
+static char report_temperature[]="\
+\\begin{figure}[!h]\n\
+\\begin{center}\n\
+\\includegraphics[width=16cm]{temperature.eps}\n\
+\\caption{Temperature over time}\n\
+\\end{center}\n\
+\\end{figure}\n\
+";
+
+static char report_end[]="\
 \\end{document}\n\
 ";
 
-static char plot_script[]="\
+static char e_plot_script[]="\
 set autoscale \n\
 unset log \n\
 unset label \n\
@@ -55,7 +76,35 @@ set xlabel 'Time [fs]' \n\
 set ylabel 'Energy per atom [eV]' \n\
 set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\
 set output 'energy.eps' \n\
-plot \"energy\" using 1:2 title 'Kinetic energy' with lines , \"energy\" using 1:3 title 'Potential energy' with lines , \"energy\" using 1:4 title 'Total energy' with lines\
+plot \"energy\" using 1:2 title 'Kinetic energy' with lines , \"energy\" using 1:3 title 'Potential energy' with lines , \"energy\" using 1:4 title 'Total energy' with lines \
+";
+
+static char pressure_plot_script[]="\
+set autoscale \n\
+unset log \n\
+unset label \n\
+set xtic auto \n\
+set ytic auto \n\
+set title 'Pressure vs. time' \n\
+set xlabel 'Time [fs]' \n\
+set ylabel 'Pressure [atm]' \n\
+set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\
+set output 'pressure.eps' \n\
+plot \"pressure\" using 1:2 title 'Pressure' with lines \
+";
+
+static char temperature_plot_script[]="\
+set autoscale \n\
+unset log \n\
+unset label \n\
+set xtic auto \n\
+set ytic auto \n\
+set title 'Temperature vs. time' \n\
+set xlabel 'Time [fs]' \n\
+set ylabel 'Temperature [K]' \n\
+set terminal postscript eps enhanced color solid lw 1 'Helvetica' 14 \n\
+set output 'temperature.eps' \n\
+plot \"temperature\" using 1:2 title 'Temperature' with lines \
 ";
 
 #endif
diff --git a/sic.c b/sic.c
index 167735f..682d89e 100644
--- a/sic.c
+++ b/sic.c
@@ -15,20 +15,51 @@
 #include "potentials/lennard_jones.h"
 #include "potentials/tersoff.h"
 
+#define INJECT		20
+#define NR_ATOMS	20	
+
 int hook(void *moldyn,void *hook_params) {
 
 	t_moldyn *md;
+	t_3dvec r,v,dist;
+	double d;
+	unsigned char run;
+	int i,j;
+	t_atom *atom;
 
 	md=moldyn;
 
- 	/* switch to direct scaling in first hook */	
-	if(md->schedule.count==0)
+	printf("\nschedule hook: ");
+
+	if(!(md->schedule.count%2)) {
+		/* add carbon at random place, and enable t scaling */
+		for(j=0;j<NR_ATOMS;j++) {
+		run=1;
+		while(run) {
+			r.x=rand_get_double(&(md->random))*md->dim.x;
+			r.y=rand_get_double(&(md->random))*md->dim.y;
+			r.z=rand_get_double(&(md->random))*md->dim.z;
+			for(i=0;i<md->count;i++) {
+				atom=&(md->atom[i]);
+				v3_sub(&dist,&(atom->r),&r);
+				d=v3_absolute_square(&dist);
+				if(d>TM_R_C)
+					run=0;
+			}
+		}
+		v.x=0; v.y=0; v.z=0;
+		add_atom(md,C,M_C,1,
+		         ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
+		         &r,&v);
+		}
+		printf("adding atoms & enable t scaling\n");
 		set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
-	/* switch off temp scaling in second hook */
-	if(md->schedule.count==1)
+	}
+	else {
+		/* disable t scaling */
+		printf("disabling t scaling\n");
 		set_pt_scale(md,0,0,0,0);
-		
-	//set_temperature(md,md->t_ref-100.0);
+	}
 
 	return 0;
 }
@@ -49,6 +80,9 @@ int main(int argc,char **argv) {
 	t_ho_params ho;
 	t_tersoff_mult_params tp;
 
+	/* atom injection counter */
+	int inject;
+
 	/* testing location & velocity vector */
 	t_3dvec r,v;
 	memset(&r,0,sizeof(t_3dvec));
@@ -176,18 +210,27 @@ int main(int argc,char **argv) {
 	thermal_init(&md,TRUE);
 
 	/* create the simulation schedule */
-	moldyn_add_schedule(&md,10001,1.0);
-	//moldyn_add_schedule(&md,501,1.0);
-	//moldyn_add_schedule(&md,501,1.0);
+	/* initial configuration */
+	moldyn_add_schedule(&md,500,1.0);
+	/* adding atoms */
+	for(inject=0;inject<INJECT;inject++) {
+		/* injecting atom and run with enabled t scaling */
+		moldyn_add_schedule(&md,400,1.0);
+		/* continue running with disabled t scaling */
+		moldyn_add_schedule(&md,100,1.0);
+	}
 
 	/* schedule hook function */
-	//moldyn_set_schedule_hook(&md,&hook,NULL);
+	moldyn_set_schedule_hook(&md,&hook,NULL);
 
 	/* activate logging */
 	moldyn_set_log_dir(&md,argv[1]);
 	moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
-	moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
+	moldyn_set_log(&md,LOG_TOTAL_ENERGY,10);
+	moldyn_set_log(&md,LOG_TEMPERATURE,10);
+	moldyn_set_log(&md,LOG_PRESSURE,10);
 	moldyn_set_log(&md,VISUAL_STEP,100);
+	moldyn_set_log(&md,SAVE_STEP,100);
 	moldyn_set_log(&md,CREATE_REPORT,0);
 
 	/*