more updates

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index f92ce55..57f76c4 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -58,6 +58,7 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
        int i;
        t_ho_params hop;
        t_lj_params ljp;
+       t_tersoff_params tp;
        double s,e;
 
        memset(moldyn,0,sizeof(t_moldyn));
@@ -81,10 +82,6 @@ int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv) {
                                        moldyn->mwrite=atoi(argv[++i]);
                                        strncpy(moldyn->mfb,argv[++i],64);
                                        break;
-                               case 'D':
-                                       moldyn->dwrite=atoi(argv[++i]);
-                                       strncpy(moldyn->dfb,argv[++i],64);
-                                       break;
                                case 'S':
                                        moldyn->swrite=atoi(argv[++i]);
                                        strncpy(moldyn->sfb,argv[++i],64);
@@ -187,16 +184,6 @@ int moldyn_log_init(t_moldyn *moldyn) {
        if(moldyn->swrite)
                moldyn->lvstat|=MOLDYN_LVSTAT_SAVE;
 
-       if(moldyn->dwrite) {
-               moldyn->dfd=open(moldyn->dfb,O_WRONLY|O_CREAT|O_TRUNC);
-                if(moldyn->dfd<0) {
-                       perror("[moldyn] dfd open");
-                       return moldyn->dfd;
-               }
-               write(moldyn->dfd,moldyn,sizeof(t_moldyn));
-               moldyn->lvstat|=MOLDYN_LVSTAT_DUMP;
-       }
-
        if((moldyn->vwrite)&&(vis)) {
                moldyn->visual=vis;
                visual_init(vis,moldyn->vfb);
@@ -544,6 +531,9 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
                }
        }
 
+       lc->dnlc=count2;
+       lc->countn=27;
+
        return count2;
 }
 
@@ -632,14 +622,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
                                        write(fd,moldyn->atom,
                                              moldyn->count*sizeof(t_atom));
                                }
+                               close(fd);
                        }       
                }
-               if(d) {
-                       if(!(i%d))
-                               write(moldyn->dfd,moldyn->atom,
-                                     moldyn->count*sizeof(t_atom));
-
-               }
                if(v) {
                        if(!(i%v)) {
                                visual_atoms(moldyn->visual,i*moldyn->tau,
@@ -687,7 +672,8 @@ printf("done\n");
 
        /* forces depending on chosen potential */
 printf("calc potential/force ...\n");
-       moldyn->potential_force_function(moldyn);
+       potential_force_calc(moldyn);
+       //moldyn->potential_force_function(moldyn);
 printf("done\n");
 
        for(i=0;i<count;i++) {
@@ -706,6 +692,100 @@ printf("done\n");
  * 
  */
 
+/* generic potential and force calculation */
+
+#define                CREATE_CELL_LIST(nb_list) \
+       link_cell_neighbour_index(moldyn,\
+                                 (atom[i].r.x+moldyn->dim.x/2)/lc->x,\
+                                 (atom[i].r.y+moldyn->dim.y/2)/lc->y,\
+                                 (atom[i].r.z+moldyn->dim.z/2)/lc->z,\
+                                 nb_list);
+
+
+int potential_force_calc(t_moldyn *moldyn) {
+
+       int i,count;
+       t_atom *atom;
+       t_linkcell *lc;
+       t_list neighbour[27];
+       t_list *this;
+       double u;
+
+       count=moldyn->count;
+       atom=moldyn->atom;
+       lc=&(moldyn->lc);
+
+       /* reset energy */
+       u=0.0;
+
+       for(i=0;i<count;i++) {
+       
+               /* reset force */
+               v3_zero(&(atom[i].f));
+
+               /* single particle potential/force */
+               if(moldyn->status&MOLDYN_STAT_1BP)
+                       moldyn->pf_func1b(moldyn,&(atom[i]));
+
+               /* 2 body pair potential/force */
+               if(moldyn->status&MOLDYN_STAT_2BP) {
+               
+                       CREATE_CELL_LIST(neighbour);
+
+                       /*
+                        * processing cell of atom i
+                        * => no need to check for empty list
+                        *    (1 element at minimum)
+                        */
+
+                       this=&(neighbour[0]);
+                       list_reset(this);
+                       do {
+                               btom=this->current->data;
+                               if(btom!=&(atom[i]))
+                                       moldyn->pf_func2b(moldyn,
+                                                         &(atom[i]),btom);
+                       } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+
+                       /*
+                        * direct neighbour cells
+                        * => no boundary condition check necessary
+                        */
+                       for(j=0;j<lc->dnlc;j++) {
+                               this=&(neighbour[j]);
+                               list_reset(this);
+                               if(this->start!=NULL) {
+                                       do {
+                                               btom=this->current->data;
+                                               moldyn->pf_func2b(moldyn,
+                                                                 &(atom[i]),
+                                                                 btom);
+                                       } while(list_next(this)!=\
+                                               L_NO_NEXT_ELEMENT);
+                       }
+
+                       /*
+                        * neighbour cells due to periodic bc
+                        * => check boundary conditions
+                        */
+                       for(j=lc->dnlc;j<lc->countn;j++) {
+                               this=&(neighbour[j]);
+                               list_reset(this);
+                               if(this->start!=NULL) {
+                                       do {
+                                               btom=this->current->data;
+                                               moldyn->pf_func2b(moldyn,
+                                                                 &(atom[i]),
+                                                                 btom);
+
+                       }
+
+               }
+
+       return 0;
+}
+
+
 /* harmonic oscillator potential and force */
 
 int harmonic_oscillator(t_moldyn *moldyn) {