added postproc stuff, sic mods (today i was a lazy ass!)

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 2fdf4c2..a310e91 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -24,6 +24,9 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 
        memset(moldyn,0,sizeof(t_moldyn));
 
+       moldyn->argc=argc;
+       moldyn->args=argv;
+
        rand_init(&(moldyn->random),NULL,1);
        moldyn->random.status|=RAND_STAT_VERBOSE;
 
@@ -510,6 +513,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
                atom[ret].brand=brand;
                atom[ret].tag=count+ret;
                check_per_bound(moldyn,&(atom[ret].r));
+               atom[ret].r_0=atom[ret].r;
        }
 
        /* update total system mass */
@@ -518,6 +522,40 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        return ret;
 }
 
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
+             t_3dvec *r,t_3dvec *v) {
+
+       t_atom *atom;
+       void *ptr;
+       int count;
+       
+       atom=moldyn->atom;
+       count=(moldyn->count)++;
+
+       ptr=realloc(atom,(count+1)*sizeof(t_atom));
+       if(!ptr) {
+               perror("[moldyn] realloc (add atom)");
+               return -1;
+       }
+       moldyn->atom=ptr;
+
+       atom=moldyn->atom;
+       atom[count].r=*r;
+       atom[count].v=*v;
+       atom[count].element=element;
+       atom[count].mass=mass;
+       atom[count].brand=brand;
+       atom[count].tag=count;
+       atom[count].attr=attr;
+       check_per_bound(moldyn,&(atom[count].r));
+       atom[count].r_0=atom[count].r;
+
+       /* update total system mass */
+       total_mass_calc(moldyn);
+
+       return 0;
+}
+
 /* cubic init */
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
@@ -630,38 +668,6 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
        return count;
 }
 
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
-             t_3dvec *r,t_3dvec *v) {
-
-       t_atom *atom;
-       void *ptr;
-       int count;
-       
-       atom=moldyn->atom;
-       count=(moldyn->count)++;
-
-       ptr=realloc(atom,(count+1)*sizeof(t_atom));
-       if(!ptr) {
-               perror("[moldyn] realloc (add atom)");
-               return -1;
-       }
-       moldyn->atom=ptr;
-
-       atom=moldyn->atom;
-       atom[count].r=*r;
-       atom[count].v=*v;
-       atom[count].element=element;
-       atom[count].mass=mass;
-       atom[count].brand=brand;
-       atom[count].tag=count;
-       atom[count].attr=attr;
-
-       /* update total system mass */
-       total_mass_calc(moldyn);
-
-       return 0;
-}
-
 int destroy_atoms(t_moldyn *moldyn) {
 
        if(moldyn->atom) free(moldyn->atom);
@@ -1097,8 +1103,10 @@ double e_kin_calc(t_moldyn *moldyn) {
        atom=moldyn->atom;
        moldyn->ekin=0.0;
 
-       for(i=0;i<moldyn->count;i++)
-               moldyn->ekin+=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+       for(i=0;i<moldyn->count;i++) {
+               atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
+               moldyn->ekin+=atom[i].ekin;
+       }
 
        return moldyn->ekin;
 }
@@ -1478,7 +1486,7 @@ return 0;
        printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
               sched->count,i,
               moldyn->t,moldyn->t_avg,
-              moldyn->p_avg/BAR,moldyn->p/BAR,
+              moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
               moldyn->volume);
        fflush(stdout);
                }
@@ -1490,13 +1498,11 @@ return 0;
        }
 
                /* check for hooks */
-               if(sched->count+1<sched->total_sched) {
-                       if(sched->hook) {
-                               printf("\n ## schedule hook %d/%d start ##\n",
-                                      sched->count+1,sched->total_sched-1);
-                               sched->hook(moldyn,sched->hook_params);
-                               printf(" ## schedule hook end ##\n");
-                       }
+               if(sched->hook) {
+                       printf("\n ## schedule hook %d/%d start ##\n",
+                              sched->count+1,sched->total_sched-1);
+                       sched->hook(moldyn,sched->hook_params);
+                       printf(" ## schedule hook end ##\n");
                }
 
                /* increase the schedule counter */