incorporated visual functions to moldyn api
[physik/posic.git] / moldyn.c
index a310e91..3c4d873 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
 #include <sys/stat.h>
 #include <fcntl.h>
 #include <unistd.h>
+#include <sys/time.h>
+#include <time.h>
 #include <math.h>
 
 #include "moldyn.h"
 #include "report/report.h"
 
+/*
+ * global variables, pse and atom colors (only needed here)
+ */ 
+
+static char *pse_name[]={
+       "*",
+       "H",
+       "He",
+       "Li",
+       "Be",
+       "B",
+       "C",
+       "N",
+       "O",
+       "F",
+       "Ne",
+       "Na",
+       "Mg",
+       "Al",
+       "Si",
+       "P",
+       "S",
+       "Cl",
+       "Ar",
+};
+
+static char *pse_col[]={
+       "*",
+       "White",
+       "He",
+       "Li",
+       "Be",
+       "B",
+       "Gray",
+       "N",
+       "Blue",
+       "F",
+       "Ne",
+       "Na",
+       "Mg",
+       "Al",
+       "Yellow",
+       "P",
+       "S",
+       "Cl",
+       "Ar",
+};
+
+/*
+ * the moldyn functions
+ */
+
 int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
 
        printf("[moldyn] init\n");
@@ -425,7 +479,6 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
                         moldyn->vlsdir);
                system(sc);
        }
-       if(&(moldyn->vis)) visual_tini(&(moldyn->vis));
 
        return 0;
 }
@@ -556,6 +609,35 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        return 0;
 }
 
+int del_atom(t_moldyn *moldyn,int tag) {
+
+       t_atom *new,*old;
+       int cnt;
+
+       old=moldyn->atom;
+
+       new=(t_atom *)malloc((moldyn->count-1)*sizeof(t_atom));
+       if(!new) {
+               perror("[moldyn]malloc (del atom)");
+               return -1;
+       }
+
+       for(cnt=0;cnt<tag;cnt++)
+               new[cnt]=old[cnt];
+       
+       for(cnt=tag+1;cnt<moldyn->count;cnt++) {
+               new[cnt-1]=old[cnt];
+               new[cnt-1].tag=cnt-1;
+       }
+
+       moldyn->count-=1;
+       moldyn->atom=new;
+
+       free(old);
+
+       return 0;
+}
+
 /* cubic init */
 int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
@@ -1348,6 +1430,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        char dir[128];
        double ds;
        double energy_scale;
+       struct timeval t1,t2;
        //double tp;
 
        sched=&(moldyn->schedule);
@@ -1368,8 +1451,8 @@ 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;
+       /* get current time */
+       gettimeofday(&t1,NULL);
 
        /* calculate initial forces */
        potential_force_calc(moldyn);
@@ -1407,6 +1490,9 @@ return 0;
                moldyn->tau_square=moldyn->tau*moldyn->tau;
                moldyn->time_steps=sched->runs[sched->count];
 
+               /* energy scaling factor (might change!) */
+               energy_scale=moldyn->count*EV;
+
        /* integration according to schedule */
 
        for(i=0;i<moldyn->time_steps;i++) {
@@ -1429,7 +1515,7 @@ return 0;
 
                /* check for log & visualization */
                if(e) {
-                       if(!(i%e))
+                       if(!(moldyn->total_steps%e))
                                dprintf(moldyn->efd,
                                        "%f %f %f %f\n",
                                        moldyn->time,moldyn->ekin/energy_scale,
@@ -1437,7 +1523,7 @@ return 0;
                                        get_total_energy(moldyn)/energy_scale);
                }
                if(m) {
-                       if(!(i%m)) {
+                       if(!(moldyn->total_steps%m)) {
                                momentum=get_total_p(moldyn);
                                dprintf(moldyn->mfd,
                                        "%f %f %f %f %f\n",moldyn->time,
@@ -1446,7 +1532,7 @@ return 0;
                        }
                }
                if(p) {
-                       if(!(i%p)) {
+                       if(!(moldyn->total_steps%p)) {
                                dprintf(moldyn->pfd,
                                        "%f %f %f %f %f\n",moldyn->time,
                                         moldyn->p/BAR,moldyn->p_avg/BAR,
@@ -1454,17 +1540,18 @@ return 0;
                        }
                }
                if(t) {
-                       if(!(i%t)) {
+                       if(!(moldyn->total_steps%t)) {
                                dprintf(moldyn->tfd,
                                        "%f %f %f\n",
                                        moldyn->time,moldyn->t,moldyn->t_avg);
                        }
                }
                if(s) {
-                       if(!(i%s)) {
+                       if(!(moldyn->total_steps%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));
@@ -1475,20 +1562,27 @@ return 0;
                        }       
                }
                if(v) {
-                       if(!(i%v)) {
+                       if(!(moldyn->total_steps%v)) {
                                visual_atoms(&(moldyn->vis),moldyn->time,
                                             moldyn->atom,moldyn->count);
                        }
                }
 
                /* display progress */
-               if(!(i%10)) {
-       printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f",
+               if(!(moldyn->total_steps%10)) {
+                       /* 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->gp_avg/BAR,
-              moldyn->volume);
+              moldyn->volume,
+              (int)(t2.tv_sec-t1.tv_sec));
        fflush(stdout);
+
+                       /* copy over time */
+                       t1=t2;
                }
 
                /* increase absolute time */
@@ -1923,6 +2017,17 @@ int moldyn_bc_check(t_moldyn *moldyn) {
        return 0;
 }
 
+/*
+ * restore function
+ */
+
+int moldyn_load(t_moldyn *moldyn) {
+
+       // later ...
+
+       return 0;
+}
+
 /*
  * post processing functions
  */
@@ -1945,3 +2050,101 @@ int get_line(int fd,char *line,int max) {
        }
 }
 
+int analyze_bonds(t_moldyn *moldyn) {
+
+       
+       
+
+       return 0;
+}
+
+/*
+ * visualization code
+ */
+
+int visual_init(t_visual *v,char *filebase) {
+
+       char file[128+8];
+
+       strncpy(v->fb,filebase,128);
+       memset(file,0,128+8);
+
+       return 0;
+}
+
+int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
+
+       int i,fd;
+       char file[128+64];
+       t_3dvec dim;
+       double help;
+
+       dim.x=v->dim.x;
+       dim.y=v->dim.y;
+       dim.z=v->dim.z;
+
+       help=(dim.x+dim.y);
+
+       sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,time);
+       fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
+       if(fd<0) {
+               perror("open visual save file fd");
+               return -1;
+       }
+
+       /* write the actual data file */
+       dprintf(fd,"# [P] %d %07.f <%f,%f,%f>\n",
+               n,time,help/40.0,help/40.0,-0.8*help);
+       for(i=0;i<n;i++)
+               dprintf(fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
+                                                atom[i].r.x,
+                                                atom[i].r.y,
+                                                atom[i].r.z,
+                                                pse_col[atom[i].element],
+                                                atom[i].ekin);
+       if(dim.x) {
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,-dim.y/2,-dim.z/2,
+                       dim.x/2,-dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,-dim.y/2,-dim.z/2,
+                       -dim.x/2,dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       dim.x/2,dim.y/2,-dim.z/2,
+                       dim.x/2,-dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,dim.y/2,-dim.z/2,
+                       dim.x/2,dim.y/2,-dim.z/2);
+
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,-dim.y/2,dim.z/2,
+                       dim.x/2,-dim.y/2,dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,-dim.y/2,dim.z/2,
+                       -dim.x/2,dim.y/2,dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       dim.x/2,dim.y/2,dim.z/2,
+                       dim.x/2,-dim.y/2,dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,dim.y/2,dim.z/2,
+                       dim.x/2,dim.y/2,dim.z/2);
+
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,-dim.y/2,dim.z/2,
+                       -dim.x/2,-dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       -dim.x/2,dim.y/2,dim.z/2,
+                       -dim.x/2,dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       dim.x/2,-dim.y/2,dim.z/2,
+                       dim.x/2,-dim.y/2,-dim.z/2);
+               dprintf(fd,"# [D] %f %f %f %f %f %f\n",
+                       dim.x/2,dim.y/2,dim.z/2,
+                       dim.x/2,dim.y/2,-dim.z/2);
+       }
+
+       close(fd);
+
+       return 0;
+}
+