added pthreads (for io only by now) + refresh screen for i count +

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index b73cd34..4788a81 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -23,6 +23,10 @@
 #include <omp.h>
 #endif
 
+#ifdef PTHREADS
+#include <pthread.h>
+#endif
+
 #include "moldyn.h"
 #include "report/report.h"
 
@@ -600,6 +604,15 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        }
        moldyn->atom=ptr;
 
+#ifdef LOWMEM_LISTS
+       ptr=realloc(moldyn->lc.subcell->list,(count+1)*sizeof(int));
+       if(!ptr) {
+               perror("[moldyn] list realloc (add atom)");
+               return -1;
+       }
+       moldyn->lc.subcell->list=ptr;
+#endif
+
        atom=moldyn->atom;
 
        /* initialize new atom */
@@ -1293,7 +1306,9 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 int link_cell_init(t_moldyn *moldyn,u8 vol) {
 
        t_linkcell *lc;
+#ifndef LOWMEM_LISTS
        int i;
+#endif
 
        lc=&(moldyn->lc);
 
@@ -1308,6 +1323,8 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 
 #ifdef STATIC_LISTS
        lc->subcell=malloc(lc->cells*sizeof(int*));
+#elif LOWMEM_LISTS
+       lc->subcell=malloc(sizeof(t_lowmem_list));
 #else
        lc->subcell=malloc(lc->cells*sizeof(t_list));
 #endif
@@ -1325,6 +1342,9 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 #ifdef STATIC_LISTS
                printf("[moldyn] initializing 'static' linked cells (%d)\n",
                       lc->cells);
+#elif LOWMEM_LISTS
+               printf("[moldyn] initializing 'lowmem' linked cells (%d)\n",
+                      lc->cells);
 #else
                printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
                       lc->cells);
@@ -1348,6 +1368,17 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
                               i,lc->subcell[0],lc->subcell);
                */
        }
+#elif LOWMEM_LISTS
+       lc->subcell->head=malloc(lc->cells*sizeof(int));
+       if(lc->subcell->head==NULL) {
+               perror("[moldyn] head init (malloc)");
+               return -1;
+       }
+       lc->subcell->list=malloc(moldyn->count*sizeof(int));
+       if(lc->subcell->list==NULL) {
+               perror("[moldyn] list init (malloc)");
+               return -1;
+       }
 #else
        for(i=0;i<lc->cells;i++)
                list_init_f(&(lc->subcell[i]));
@@ -1362,22 +1393,26 @@ int link_cell_init(t_moldyn *moldyn,u8 vol) {
 int link_cell_update(t_moldyn *moldyn) {
 
        int count,i,j,k;
-       int nx,ny;
+       int nx,nxy;
        t_atom *atom;
        t_linkcell *lc;
 #ifdef STATIC_LISTS
        int p;
+#elif LOWMEM_LISTS
+       int p;
 #endif
 
        atom=moldyn->atom;
        lc=&(moldyn->lc);
 
        nx=lc->nx;
-       ny=lc->ny;
+       nxy=nx*lc->ny;
 
        for(i=0;i<lc->cells;i++)
 #ifdef STATIC_LISTS
-               memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+               memset(lc->subcell[i],-1,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#elif LOWMEM_LISTS
+               lc->subcell->head[i]=-1;
 #else
                list_destroy_f(&(lc->subcell[i]));
 #endif
@@ -1389,7 +1424,7 @@ int link_cell_update(t_moldyn *moldyn) {
        
 #ifdef STATIC_LISTS
                p=0;
-               while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
+               while(lc->subcell[i+j*nx+k*nxy][p]!=-1)
                        p++;
 
                if(p>=MAX_ATOMS_PER_LIST) {
@@ -1397,9 +1432,13 @@ int link_cell_update(t_moldyn *moldyn) {
                        return -1;
                }
 
-               lc->subcell[i+j*nx+k*nx*ny][p]=count;
+               lc->subcell[i+j*nx+k*nxy][p]=count;
+#elif LOWMEM_LISTS
+               p=i+j*nx+k*nxy;
+               lc->subcell->list[count]=lc->subcell->head[p];
+               lc->subcell->head[p]=count;
 #else
-               list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
+               list_add_immediate_f(&(lc->subcell[i+j*nx+k*nxy]),
                                     &(atom[count]));
                /*
                if(j==0&&k==0)
@@ -1415,6 +1454,8 @@ int link_cell_update(t_moldyn *moldyn) {
 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 #ifdef STATIC_LISTS
                               int **cell
+#elif LOWMEM_LISTS
+                              int *cell
 #else
                               t_list *cell
 #endif
@@ -1440,7 +1481,11 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
                printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
                       i,nx,j,ny,k,nz);
 
+#ifndef LOWMEM_LISTS
        cell[0]=lc->subcell[i+j*nx+k*a];
+#else
+       cell[0]=lc->subcell->head[i+j*nx+k*a];
+#endif
        for(ci=-1;ci<=1;ci++) {
                bx=0;
                x=i+ci;
@@ -1464,10 +1509,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
                                }
                                if(!(ci|cj|ck)) continue;
                                if(bx|by|bz) {
+#ifndef LOWMEM_LISTS
                                        cell[--count2]=lc->subcell[x+y*nx+z*a];
+#else
+                               cell[--count2]=lc->subcell->head[x+y*nx+z*a];
+#endif
+                                       
                                }
                                else {
+#ifndef LOWMEM_LISTS
                                        cell[count1++]=lc->subcell[x+y*nx+z*a];
+#else
+                               cell[count1++]=lc->subcell->head[x+y*nx+z*a];
+#endif
                                }
                        }
                }
@@ -1480,11 +1534,19 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
 
 int link_cell_shutdown(t_moldyn *moldyn) {
 
+#ifndef LOWMEM_LISTS
        int i;
+#endif
        t_linkcell *lc;
 
        lc=&(moldyn->lc);
 
+#if LOWMEM_LISTS
+       free(lc->subcell->head);
+       free(lc->subcell->list);
+
+#else
+
        for(i=0;i<lc->cells;i++) {
 #ifdef STATIC_LISTS
                free(lc->subcell[i]);
@@ -1493,6 +1555,7 @@ int link_cell_shutdown(t_moldyn *moldyn) {
                list_destroy_f(&(lc->subcell[i]));
 #endif
        }
+#endif
 
        free(lc->subcell);
 
@@ -1545,6 +1608,33 @@ int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params) {
  *
  */
 
+/* helper / special functions */
+
+#ifdef PTHREADS
+void *write_save_file(void *ptr) {
+
+       int fd;
+       char dir[128];
+       t_moldyn *moldyn;
+
+       moldyn=ptr;
+
+       snprintf(dir,128,"%s/s-%08.f.save",moldyn->vlsdir,moldyn->time);
+       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));
+               write(fd,moldyn->atom,moldyn->count*sizeof(t_atom));
+       }
+
+       close(fd);
+
+       pthread_exit(NULL);
+
+       return 0;
+}
+#endif
+
 /* start the integration */
 
 int moldyn_integrate(t_moldyn *moldyn) {
@@ -1554,13 +1644,26 @@ int moldyn_integrate(t_moldyn *moldyn) {
        t_3dvec momentum;
        t_moldyn_schedule *sched;
        t_atom *atom;
+#ifndef PTHREADS
        int fd;
        char dir[128];
+#endif
        double ds;
        double energy_scale;
        struct timeval t1,t2;
        //double tp;
 
+#ifdef PTHREADS
+       u8 first,change;
+       pthread_t io_thread;
+       int ret;
+       t_moldyn md_copy;
+       t_atom *atom_copy;
+
+       first=1;
+       change=0;
+#endif
+
        sched=&(moldyn->schedule);
        atom=moldyn->atom;
 
@@ -1609,6 +1712,16 @@ int moldyn_integrate(t_moldyn *moldyn) {
        /* debugging, ignore */
        moldyn->debug=0;
 
+       /* zero & init moldyn copy */
+#ifdef PTHREADS
+       memset(&md_copy,0,sizeof(t_moldyn));
+       atom_copy=malloc(moldyn->count*sizeof(t_atom));
+       if(atom_copy==NULL) {
+               perror("[moldyn] malloc atom copy (init)");
+               return -1;
+       }
+#endif
+
        /* tell the world */
        printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1693,7 +1806,35 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
                if(s) {
                        if(!(moldyn->total_steps%s)) {
-                               snprintf(dir,128,"%s/s-%07.f.save",
+#ifdef PTHREADS
+       /* check whether thread has not terminated yet */
+       if(!first) {
+               ret=pthread_join(io_thread,NULL);
+       }
+       first=0;
+       /* prepare and start thread */
+       if(moldyn->count!=md_copy.count) {
+               free(atom_copy);
+               change=1;
+       }
+       memcpy(&md_copy,moldyn,sizeof(t_moldyn));
+       if(change) {
+               atom_copy=malloc(moldyn->count*sizeof(t_atom));
+               if(atom_copy==NULL) {
+                       perror("[moldyn] malloc atom copy (change)");
+                       return -1;
+               }
+       }
+       md_copy.atom=atom_copy;
+       memcpy(atom_copy,moldyn->atom,moldyn->count*sizeof(t_atom));
+       change=0;
+       ret=pthread_create(&io_thread,NULL,write_save_file,&md_copy);
+       if(ret) {
+               perror("[moldyn] create save file thread\n");
+               return -1;
+       }
+#else
+                               snprintf(dir,128,"%s/s-%08.f.save",
                                         moldyn->vlsdir,moldyn->time);
                                fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT,
                                        S_IRUSR|S_IWUSR);
@@ -1704,6 +1845,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                                              moldyn->count*sizeof(t_atom));
                                }
                                close(fd);
+#endif
                        }       
                }
                if(a) {
@@ -1713,7 +1855,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
 
                /* display progress */
-               //if(!(moldyn->total_steps%10)) {
+               if(!(i%10)) {
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
@@ -1729,7 +1871,7 @@ printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
 
                        /* copy over time */
                        t1=t2;
-               //}
+               }
 
                /* increase absolute time */
                moldyn->time+=moldyn->tau;
@@ -1828,6 +1970,9 @@ int potential_force_calc(t_moldyn *moldyn) {
        int *neighbour_i[27];
        int p,q;
        t_atom *atom;
+#elif LOWMEM_LISTS
+       int neighbour_i[27];
+       int p,q;
 #else
        t_list neighbour_i[27];
        t_list neighbour_i2[27];
@@ -1851,6 +1996,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
        /* reset force, site energy and virial of every atom */
 #ifdef PARALLEL
+       i=omp_get_thread_num();
        #pragma omp parallel for private(virial)
 #endif
        for(i=0;i<count;i++) {
@@ -1903,10 +2049,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 #ifdef STATIC_LISTS
                                p=0;
 
-                               while(neighbour_i[j][p]!=0) {
+                               while(neighbour_i[j][p]!=-1) {
 
                                        jtom=&(atom[neighbour_i[j][p]]);
                                        p++;
+#elif LOWMEM_LISTS
+                               p=neighbour_i[j];
+
+                               while(p!=-1) {
+
+                                       jtom=&(itom[p]);
+                                       p=lc->subcell->list[p];
 #else
                                this=&(neighbour_i[j]);
                                list_reset_f(this);
@@ -1930,6 +2083,8 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        }
 #ifdef STATIC_LISTS
                                }
+#elif LOWMEM_LISTS
+                               }
 #else
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
@@ -1945,6 +2100,8 @@ int potential_force_calc(t_moldyn *moldyn) {
                /* copy the neighbour lists */
 #ifdef STATIC_LISTS
                /* no copy needed for static lists */
+#elif LOWMEM_LISTS
+               /* no copy needed for lowmem lists */
 #else
                memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
 #endif
@@ -1956,10 +2113,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 #ifdef STATIC_LISTS
                        p=0;
 
-                       while(neighbour_i[j][p]!=0) {
+                       while(neighbour_i[j][p]!=-1) {
 
                                jtom=&(atom[neighbour_i[j][p]]);
                                p++;
+#elif LOWMEM_LISTS
+                               p=neighbour_i[j];
+
+                               while(p!=-1) {
+
+                                       jtom=&(itom[p]);
+                                       p=lc->subcell->list[p];
 #else
                        this=&(neighbour_i[j]);
                        list_reset_f(this);
@@ -2000,10 +2164,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 #ifdef STATIC_LISTS
                                        q=0;
 
-                                       while(neighbour_i[j][q]!=0) {
+                                       while(neighbour_i[k][q]!=-1) {
 
                                                ktom=&(atom[neighbour_i[k][q]]);
                                                q++;
+#elif LOWMEM_LISTS
+                                       q=neighbour_i[k];
+
+                                       while(q!=-1) {
+
+                                               ktom=&(itom[q]);
+                                               q=lc->subcell->list[q];
 #else
                                        that=&(neighbour_i2[k]);
                                        list_reset_f(that);
@@ -2031,6 +2202,8 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                                  bc_ik|bc_ij);
 #ifdef STATIC_LISTS
                                        }
+#elif LOWMEM_LISTS
+                                       }
 #else
                                        } while(list_next_f(that)!=\
                                                L_NO_NEXT_ELEMENT);
@@ -2055,10 +2228,17 @@ int potential_force_calc(t_moldyn *moldyn) {
 #ifdef STATIC_LISTS
                                        q=0;
 
-                                       while(neighbour_i[j][q]!=0) {
+                                       while(neighbour_i[k][q]!=-1) {
 
                                                ktom=&(atom[neighbour_i[k][q]]);
                                                q++;
+#elif LOWMEM_LISTS
+                                       q=neighbour_i[k];
+
+                                       while(q!=-1) {
+
+                                               ktom=&(itom[q]);
+                                               q=lc->subcell->list[q];
 #else
                                        that=&(neighbour_i2[k]);
                                        list_reset_f(that);
@@ -2087,6 +2267,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
 #ifdef STATIC_LISTS
                                        }
+#elif LOWMEM_LISTS
+                                       }
 #else
                                        } while(list_next_f(that)!=\
                                                L_NO_NEXT_ELEMENT);
@@ -2104,6 +2286,8 @@ int potential_force_calc(t_moldyn *moldyn) {
                                }
 #ifdef STATIC_LISTS
                        }
+#elif LOWMEM_LISTS
+                       }
 #else
                        } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
@@ -2345,6 +2529,9 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 #ifdef STATIC_LISTS
        int *neighbour[27];
        int p;
+#elif LOWMEM_LISTS
+       int neighbour[27];
+       int p;
 #else
        t_list neighbour[27];
        t_list *this;
@@ -2371,10 +2558,17 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 #ifdef STATIC_LISTS
                        p=0;
 
-                       while(neighbour[j][p]!=0) {
+                       while(neighbour[j][p]!=-1) {
 
                                jtom=&(moldyn->atom[neighbour[j][p]]);
                                p++;
+#elif LOWMEM_LISTS
+                       p=neighbour[j];
+
+                       while(p!=-1) {
+
+                               jtom=&(itom[p]);
+                               p=lc->subcell->list[p];
 #else
                        this=&(neighbour[j]);
                        list_reset_f(this);
@@ -2392,6 +2586,8 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 #ifdef STATIC_LISTS
                        }
+#elif LOWMEM_LISTS
+                       }
 #else
                        } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
@@ -2402,6 +2598,84 @@ int process_2b_bonds(t_moldyn *moldyn,void *data,
 
 }
 
+/*
+ * function to find neighboured atoms
+ */
+
+int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom,
+                       int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom,
+                                     void *data,u8 bc)) {
+
+       t_linkcell *lc;
+#ifdef STATIC_LISTS
+       int *neighbour[27];
+       int p;
+#elif LOWMEM_LISTS
+       int neighbour[27];
+       int p;
+#else
+       t_list neighbour[27];
+       t_list *this;
+#endif
+       u8 bc;
+       t_atom *natom;
+       int j;
+
+       lc=&(moldyn->lc);
+       
+       /* neighbour indexing */
+       link_cell_neighbour_index(moldyn,
+                                 (atom->r.x+moldyn->dim.x/2)/lc->x,
+                                 (atom->r.y+moldyn->dim.y/2)/lc->x,
+                                 (atom->r.z+moldyn->dim.z/2)/lc->x,
+                                 neighbour);
+
+       for(j=0;j<27;j++) {
+
+               bc=(j<lc->dnlc)?0:1;
+
+#ifdef STATIC_LISTS
+               p=0;
+
+               while(neighbour[j][p]!=-1) {
+
+                       natom=&(moldyn->atom[neighbour[j][p]]);
+                       p++;
+#elif LOWMEM_LISTS
+               p=neighbour[j];
+
+               while(p!=-1) {
+
+                       natom=&(moldyn->atom[p]);
+                       p=lc->subcell->list[p];
+#else
+               this=&(neighbour[j]);
+               list_reset_f(this);
+
+               if(this->start==NULL)
+                       continue;
+
+               do {
+
+                       natom=this->current->data;
+#endif
+
+                       /* process bond */
+                       process(moldyn,atom,natom,data,bc);
+
+#ifdef STATIC_LISTS
+               }
+#elif LOWMEM_LISTS
+               }
+#else
+               } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
+#endif
+       }
+
+       return 0;
+
+}
+
 /*
  * post processing functions
  */
@@ -2746,7 +3020,7 @@ int visual_atoms(t_moldyn *moldyn) {
 
        help=(dim.x+dim.y);
 
-       sprintf(file,"%s/atomic_conf_%07.f.xyz",v->fb,moldyn->time);
+       sprintf(file,"%s/atomic_conf_%08.f.xyz",v->fb,moldyn->time);
        vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
        if(vb.fd<0) {
                perror("open visual save file fd");
@@ -2756,7 +3030,7 @@ int visual_atoms(t_moldyn *moldyn) {
        /* write the actual data file */
 
        // povray header
-       dprintf(vb.fd,"# [P] %d %07.f <%f,%f,%f>\n",
+       dprintf(vb.fd,"# [P] %d %08.f <%f,%f,%f>\n",
                moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
 
        // atomic configuration