Makefile: clean everything, moldyn: add time info

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 992d112..92c4e58 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -13,6 +13,8 @@
 #include <sys/stat.h>
 #include <fcntl.h>
 #include <unistd.h>
+#include <sys/time.h>
+#include <time.h>
 #include <math.h>
 
 #include "moldyn.h"
@@ -24,6 +26,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 +515,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 +524,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 +670,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 +1105,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;
 }
@@ -1340,6 +1350,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        char dir[128];
        double ds;
        double energy_scale;
+       struct timeval t1,t2;
        //double tp;
 
        sched=&(moldyn->schedule);
@@ -1363,6 +1374,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* energy scaling factor */
        energy_scale=moldyn->count*EV;
 
+       /* get current time */
+       gettimeofday(&t1,NULL);
+
        /* calculate initial forces */
        potential_force_calc(moldyn);
 #ifdef DEBUG
@@ -1391,7 +1405,8 @@ return 0;
        printf("[moldyn] integration start, go get a coffee ...\n");
 
        /* executing the schedule */
-       for(sched->count=0;sched->count<sched->total_sched;sched->count++) {
+       sched->count=0;
+       while(sched->count<sched->total_sched) {
 
                /* setting amount of runs and finite time step size */
                moldyn->tau=sched->tau[sched->count];
@@ -1455,7 +1470,8 @@ return 0;
                        if(!(i%s)) {
                                snprintf(dir,128,"%s/s-%07.f.save",
                                         moldyn->vlsdir,moldyn->time);
-                               fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
+                               fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
+                                       S_IRUSR|S_IWUSR);
                                if(fd<0) perror("[moldyn] save fd open");
                                else {
                                        write(fd,moldyn,sizeof(t_moldyn));
@@ -1474,12 +1490,19 @@ return 0;
 
                /* display progress */
                if(!(i%10)) {
-       printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
+                       /* get current time */
+                       gettimeofday(&t2,NULL);
+
+       printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
               sched->count,i,
               moldyn->t,moldyn->t_avg,
-              moldyn->p_avg/BAR,moldyn->p/BAR,
-              moldyn->volume);
+              moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
+              moldyn->volume,
+              (int)(t2.tv_sec-t1.tv_sec));
        fflush(stdout);
+
+                       /* copy over time */
+                       t1=t2;
                }
 
                /* increase absolute time */
@@ -1489,13 +1512,15 @@ 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);
-                               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 */
+               sched->count+=1;
 
        }
 
@@ -1912,6 +1937,17 @@ int moldyn_bc_check(t_moldyn *moldyn) {
        return 0;
 }
 
+/*
+ * restore function
+ */
+
+int moldyn_load(t_moldyn *moldyn) {
+
+       // later ...
+
+       return 0;
+}
+
 /*
  * post processing functions
  */
@@ -1934,3 +1970,11 @@ int get_line(int fd,char *line,int max) {
        }
 }
 
+int analyze_bonds(t_moldyn *moldyn) {
+
+       
+       
+
+       return 0;
+}
+