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

testing

diff --git a/moldyn.c b/moldyn.c
index ca6f84a65e858b8f21534cb853f7b85eb81460fc..9a62b9f7c93d5142f5b15d57d099f4ba160604ce 100644 (file)
--- 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 b14037219f7ddcdf5b34adf4deb173d78636510a..85a4950cc513c9cde9b7dc0cb66f9d49503e3b8f 100644 (file)
--- 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 01364ea75aa587f37d829bd335566679dad0491d..16f771c7be4b614f10c6fd70b9476fb0bfd5e23f 100644 (file)
--- 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 e46ff2a6e55a2b4a1dc8aa668f6d9ebd8feabe14..4fe42f77819c782a38881d10c9ac348af56183c9 100644 (file)
--- 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 167735f2fa75f7ebdcc32d12b5d4f9f66c61dc05..682d89ee6fcd46d96b1cb341b6c4b8d708a61a85 100644 (file)
--- 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);
 
        /*