added report system, still messing around with virial pressure

diff --git a/clean b/clean
index 1f359db..e702e9e 100755 (executable)
--- a/clean
+++ b/clean
@@ -1,2 +1,2 @@
-rm -rf $1
+rm -rf $1/*
 mkdir -p $1

diff --git a/moldyn.c b/moldyn.c
index 94d596e..817a727 100644 (file)
--- 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;

diff --git a/moldyn.h b/moldyn.h
index c910d7f..7329099 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -12,6 +12,7 @@
 #include "random/random.h"
 #include "list/list.h"
 
+#include "report/report.h"
 
 /*
  * 
@@ -82,7 +83,6 @@ typedef struct s_moldyn {
 
        t_3dvec dim;            /* dimensions of the simulation volume */
        double volume;          /* volume of sim cell (dim.x*dim.y*dim.z) */
-       double vt1,vt2;
 
        /* potential force function and parameter pointers */
        int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
@@ -139,6 +139,10 @@ typedef struct s_moldyn {
        int mfd;                /* fd for momentum 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 */
 
        u8 status;              /* general moldyn properties */
 
@@ -162,7 +166,6 @@ typedef struct s_moldyn {
 #define P_SCALE_BERENDSEN              0x04    /* berendsen p control */
 #define P_SCALE_DIRECT                 0x08    /* direct p control */
 
-
 /*
  *
  * potential parameter structures
@@ -307,6 +310,7 @@ typedef struct s_tersoff_mult_params {
 #define LOG_TOTAL_MOMENTUM             0x02
 #define SAVE_STEP                      0x04
 #define VISUAL_STEP                    0x08
+#define CREATE_REPORT                  0x10
 
 #define TRUE                           1
 #define FALSE                          0
@@ -361,8 +365,9 @@ typedef struct s_tersoff_mult_params {
  * lattice constants
  */
 
-#define FCC                    0x01
-#define DIAMOND                        0x02
+#define CUBIC                  0x01
+#define FCC                    0x02
+#define DIAMOND                        0x04
 
 
 /*
@@ -393,11 +398,13 @@ int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params);
 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params);
 
 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
+int moldyn_set_report(t_moldyn *moldyn,char *author,char *title);
 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
 int moldyn_log_shutdown(t_moldyn *moldyn);
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                    u8 attr,u8 brand,int a,int b,int c);
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
 int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,

diff --git a/report/report.h b/report/report.h
new file mode 100644 (file)
index 0000000..7fd3180
--- /dev/null
+++ b/report/report.h
@@ -0,0 +1,64 @@
+/*
+ * report.h - header file of the report subsystem
+ *
+ * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
+ *
+ */
+
+#ifndef REPORT_H
+#define REPORT_H
+
+static char report_start[]="\
+\\pdfoutput=0\n\
+\\documentclass[a4paper,11pt]{report}\n\
+\\usepackage[activate]{pdfcprot}\n\
+\\usepackage{verbatim}\n\
+\\usepackage{a4}\n\
+\\usepackage{a4wide}\n\
+\\usepackage[english,german]{babel}\n\
+\\usepackage[latin1]{inputenc}\n\
+\\usepackage[T1]{fontenc}\n\
+\\usepackage{amsmath}\n\
+\\usepackage{ae}\n\
+\\usepackage{aecompl}\n\
+\\usepackage[dvips]{graphicx}\n\
+\\graphicspath{{./}}\n\
+\\usepackage{color}\n\
+\\usepackage{pstricks}\n\
+\\usepackage{pst-node}\n\
+\\usepackage{rotating}\n\
+\\selectlanguage{english}\n\
+\\author{%s}\n\
+\\title{%s}\n\
+\\begin{document}\n\
+\\maketitle\n\
+";
+
+static char report_end[]="\
+\\begin{figure}[!h]\n\
+\\begin{center}\n\
+\\includegraphics[width=10cm]{energy.eps}\n\
+\\caption{Kinetic, potential and total energy over time}\n\
+\\end{center}\n\
+\\end{figure}\n\
+\\end{document}\n\
+";
+
+static char plot_script[]="\
+set autoscale \n\
+unset log \n\
+unset label \n\
+set xtic auto \n\
+set ytic auto \n\
+set title 'Energy vs. time' \n\
+set xlabel 'Time [fs]' \n\
+set ylabel 'Energy [eV]' \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 \n\
+#set size 1.0, 0.6 \n\
+set terminal postscript eps enhanced color dashed lw 1 'Helvetica' 14 \n\
+set output 'energy.eps' \n\
+replot\
+";
+
+
+#endif

diff --git a/sic.c b/sic.c
index 5cfbb03..b5d4596 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -24,10 +24,8 @@ int hook(void *moldyn,void *hook_params) {
        if(md->schedule.count==1)
                set_pt_scale(md,0,0,0,0);
                
-       
        //set_temperature(md,md->t_ref-100.0);
 
-
        return 0;
 }
 
@@ -66,7 +64,7 @@ int main(int argc,char **argv) {
 
        /* cutoff radius */
        //set_cutoff(&md,TM_S_SI);
-       set_cutoff(&md,2*LC_SI);
+       set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5));
 
        /*
         * potential parameters
@@ -115,16 +113,17 @@ int main(int argc,char **argv) {
        tersoff_mult_complete_params(&tp);
 
        /* set (initial) dimensions of simulation volume */
-       set_dim(&md,6*LC_SI,6*LC_SI,6*LC_SI,TRUE);
+       set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE);
 
        /* set periodic boundary conditions in all directions */
        set_pbc(&md,TRUE,TRUE,TRUE);
 
        /* create the lattice / place atoms */
-       create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
+       create_lattice(&md,FCC,LC_SI*0.5*sqrt(1.5),SI,M_SI,
+       //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
        //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
                       ATOM_ATTR_2BP|ATOM_ATTR_HB,
-                      0,6,6,6);
+                      0,8,8,8);
        moldyn_bc_check(&md);
 
        /* testing configuration */
@@ -161,24 +160,26 @@ int main(int argc,char **argv) {
        /* set p/t scaling */
        //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
        //                 T_SCALE_BERENDSEN,100.0);
-       set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
+       //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
        //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0);
        
        /* initial thermal fluctuations of particles (in equilibrium) */
        thermal_init(&md,TRUE);
 
        /* create the simulation schedule */
-       moldyn_add_schedule(&md,1001,3.0);
-       moldyn_add_schedule(&md,501,1.0);
-       moldyn_add_schedule(&md,501,1.0);
+       moldyn_add_schedule(&md,1001,1.0);
+       //moldyn_add_schedule(&md,501,1.0);
+       //moldyn_add_schedule(&md,501,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_log(&md,LOG_TOTAL_ENERGY,5);
-       moldyn_set_log(&md,VISUAL_STEP,50);
+       moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
+       moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
+       moldyn_set_log(&md,VISUAL_STEP,10);
+       moldyn_set_log(&md,CREATE_REPORT,0);
 
        /*
         * let's do the actual md algorithm now

diff --git a/visual/visual.c b/visual/visual.c
index a6fda4f..a76c608 100644 (file)
--- a/visual/visual.c
+++ b/visual/visual.c
@@ -87,7 +87,7 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
 
        /* script file update */
        dprintf(v->fd,"load xyz %s\n",file);
-       dprintf(v->fd,"spacefill 270\n");
+       dprintf(v->fd,"spacefill 100\n");
        dprintf(v->fd,"rotate x 100\n");
        dprintf(v->fd,"rotate y 10\n");
        dprintf(v->fd,"set ambient 20\n");