security ci

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index f92ce55..d74b391 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,9 +692,125 @@ printf("done\n");
  * 
  */
 
+/* generic potential and force calculation */
+
+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;
+       unsigned char bc,bc3;
+       int countn,dnlc;
+
+       count=moldyn->count;
+       atom=moldyn->atom;
+       lc=&(moldyn->lc);
+
+       /* reset energy */
+       moldyn->energy=0.0;
+
+       for(i=0;i<count;i++) {
+       
+               /* reset force */
+               v3_zero(&(atom[i].f));
+
+               /* single particle potential/force */
+               if(atom[i].attr&ATOM_ATTR_1BP)
+                       moldyn->pf_func1b(moldyn,&(atom[i]));
+
+               /* 2 body pair potential/force */
+               if(atom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
+               
+                       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,
+                               neighbour);
+
+                       countn=lc->countn;
+                       dnlc=lc->dnlc;
+
+                       for(j=0;j<countn;j++) {
+
+                               this=&(neighbour[j]);
+                               list_reset(this);
+
+                               if(this->start==NULL)
+                                       continue;
+
+                               bc=(j<dnlc)?0:1;
+
+                               do {
+                                       btom=this->current->data;
+
+                                       if(btom==&(atom[i]))
+                                               continue;
+
+                                       if((btom->attr&ATOM_ATTR_2BP)&
+                                          (atom[i].attr&ATOM_ATTR_2BP))
+                                               moldyn->pf_func2b(moldyn,
+                                                                 &(atom[i]),
+                                                                 btom,
+                                                                 bc);
+
+                                       /* 3 body potential/force */
+
+                                       if(!(atom[i].attr&ATOM_ATTR_3BP)||
+                                          !(btom->attr&ATOM_ATTR_3BP))
+                                               continue;
+
+                                       link_cell_neighbour_index(moldyn,
+                                          (btom->r.x+moldyn->dim.x/2)/lc->x,
+                                          (btom->r.y+moldyn->dim.y/2)/lc->y,
+                                          (btom->r.z+moldyn->dim.z/2)/lc->z,
+                                          neighbourk);
+
+                                       for(k=0;k<lc->countn;k++) {
+
+                                               thisk=&(neighbourk[k]);
+                                               list_reset(thisk);
+                                       
+                                               if(thisk->start==NULL)
+                                                       continue;
+
+                                               bck=(k<lc->dnlc)?0:1;
+
+                                               do {
+
+                       ktom=thisk->current->data;
+
+                       if(!(ktom->attr&ATOM_ATTR_3BP))
+                               continue;
+
+                       if(ktom==btom)
+                               continue;
+
+                       if(ktom==&(atom[i]))
+                               continue;
+
+                       moldyn->pf_func3b(moldyn,&(atom[i]),btom,ktom,bck);
+
+                                               } while(list_next(thisk)!=\
+                                                       L_NO_NEXT_ELEMENT);
+                                       
+                               } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+                       }
+               }
+       }
+
+       return 0;
+}
+
+/*
+ * example potentials
+ */
+
 /* harmonic oscillator potential and force */
 
-int harmonic_oscillator(t_moldyn *moldyn) {
+int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,unsigned char bc)) {
 
        t_ho_params *params;
        t_atom *atom,*btom;