]> hackdaworld.org Git - physik/posic.git/commitdiff
safety checkin
authorhackbard <hackbard>
Tue, 12 Dec 2006 11:50:04 +0000 (11:50 +0000)
committerhackbard <hackbard>
Tue, 12 Dec 2006 11:50:04 +0000 (11:50 +0000)
clean
moldyn.c
moldyn.h
ppm2avi
run
sic.c
visual/visual.c

diff --git a/clean b/clean
index ff74c9b029ad16c05411ccb5d8329b608063e5b9..ca57788d1f2cde4c2146b3f22857ab92a2dcc22e 100755 (executable)
--- a/clean
+++ b/clean
@@ -1 +1 @@
-rm -f saves/* video/*.ppm video/*.png > /dev/null 2>&1
+rm -f $1/* > /dev/null 2>&1
index 9d9901c55d90c3effd806ded40ce2d93dc1b2430..7ea5d5080ab1ab185fa46b2137894239e9233902 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -155,38 +155,56 @@ int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params) {
        return 0;
 }
 
-int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer) {
+int moldyn_set_log_dir(t_moldyn *moldyn,char *dir) {
+
+       strncpy(moldyn->vlsdir,dir,127);
+
+       return 0;
+}
+       
+int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer) {
+
+       char filename[128];
+       int ret;
 
        switch(type) {
                case LOG_TOTAL_ENERGY:
                        moldyn->ewrite=timer;
-                       moldyn->efd=open(fb,O_WRONLY|O_CREAT|O_TRUNC);
+                       snprintf(filename,127,"%s/energy",moldyn->vlsdir);
+                       moldyn->efd=open(filename,
+                                        O_WRONLY|O_CREAT|O_EXCL,
+                                        S_IRUSR|S_IWUSR);
                        if(moldyn->efd<0) {
-                               perror("[moldyn] efd open");
+                               perror("[moldyn] energy log fd open");
                                return moldyn->efd;
                        }
                        dprintf(moldyn->efd,"# total energy log file\n");
                        break;
                case LOG_TOTAL_MOMENTUM:
                        moldyn->mwrite=timer;
-                       moldyn->mfd=open(fb,O_WRONLY|O_CREAT|O_TRUNC);
+                       snprintf(filename,127,"%s/momentum",moldyn->vlsdir);
+                       moldyn->mfd=open(filename,
+                                        O_WRONLY|O_CREAT|O_EXCL,
+                                        S_IRUSR|S_IWUSR);
                        if(moldyn->mfd<0) {
-                               perror("[moldyn] mfd open");
+                               perror("[moldyn] momentum log fd open");
                                return moldyn->mfd;
                        }
                        dprintf(moldyn->efd,"# total momentum log file\n");
                        break;
                case SAVE_STEP:
                        moldyn->swrite=timer;
-                       strncpy(moldyn->sfb,fb,63);
                        break;
                case VISUAL_STEP:
                        moldyn->vwrite=timer;
-                       strncpy(moldyn->vfb,fb,63);
-                       visual_init(&(moldyn->vis),fb);
+                       ret=visual_init(&(moldyn->vis),moldyn->vlsdir);
+                       if(ret<0) {
+                               printf("[moldyn] visual init failure\n");
+                               return ret;
+                       }
                        break;
                default:
-                       printf("unknown log mechanism: %02x\n",type);
+                       printf("[moldyn] unknown log mechanism: %02x\n",type);
                        return -1;
        }
 
@@ -371,8 +389,10 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
        /* (temporary) hack for e,t = 0 */
        if(e==0.0) {
        moldyn->t=0.0;
-               if(moldyn->t_ref!=0.0)
+               if(moldyn->t_ref!=0.0) {
                        thermal_init(moldyn,equi_init);
+                       return 0;
+               }
                else
                        return 0; /* no scaling needed */
        }
@@ -384,17 +404,13 @@ int scale_velocity(t_moldyn *moldyn,u8 equi_init) {
                scale*=2.0;
        else
                if(moldyn->pt_scale&T_SCALE_BERENDSEN)
-                       scale=1.0+moldyn->tau*(scale-1.0)/moldyn->t_tc;
-printf("scale=%f\n",scale);
+                       scale=1.0+(scale-1.0)/moldyn->t_tc;
        scale=sqrt(scale);
-printf("debug: %f %f %f %f \n",scale,moldyn->t_ref,moldyn->t,moldyn->tau);
 
        /* velocity scaling */
        for(i=0;i<moldyn->count;i++) {
-printf("vorher: %f\n",atom[i].v.x);
                if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
                        v3_scale(&(atom[i].v),&(atom[i].v),scale);
-printf("nachher: %f\n",atom[i].v.x);
        }
 
        return 0;
@@ -644,7 +660,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        t_moldyn_schedule *schedule;
        t_atom *atom;
        int fd;
-       char fb[128];
+       char dir[128];
        double ds;
 
        schedule=&(moldyn->schedule);
@@ -702,14 +718,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                if(moldyn->pt_scale&(T_SCALE_BERENDSEN|T_SCALE_DIRECT))
                        scale_velocity(moldyn,FALSE);
 
-               /* increase absolute time */
-               moldyn->time+=moldyn->tau;
-
                /* check for log & visualization */
                if(e) {
                        if(!(i%e))
                                dprintf(moldyn->efd,
-                                       "%.15f %.45f %.45f %.45f\n",
+                                       "%f %f %f %f\n",
                                        moldyn->time,update_e_kin(moldyn),
                                        moldyn->energy,
                                        get_total_energy(moldyn));
@@ -718,15 +731,14 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        if(!(i%m)) {
                                p=get_total_p(moldyn);
                                dprintf(moldyn->mfd,
-                                       "%.15f %.45f\n",moldyn->time,
-                                       v3_norm(&p));
+                                       "%f %f\n",moldyn->time,v3_norm(&p));
                        }
                }
                if(s) {
                        if(!(i%s)) {
-                               snprintf(fb,128,"%s-%f-%.15f.save",moldyn->sfb,
-                                        moldyn->t,i*moldyn->tau);
-                               fd=open(fb,O_WRONLY|O_TRUNC|O_CREAT);
+                               snprintf(dir,128,"%s/s-%07.f.save",
+                                        moldyn->vlsdir,moldyn->time);
+                               fd=open(dir,O_WRONLY|O_TRUNC|O_CREAT);
                                if(fd<0) perror("[moldyn] save fd open");
                                else {
                                        write(fd,moldyn,sizeof(t_moldyn));
@@ -746,6 +758,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        }
                }
 
+               /* increase absolute time */
+               moldyn->time+=moldyn->tau;
+
        }
 
                /* check for hooks */
@@ -812,7 +827,6 @@ int velocity_verlet(t_moldyn *moldyn) {
 /* generic potential and force calculation */
 
 int potential_force_calc(t_moldyn *moldyn) {
-printf("start pot force calc\n");
 
        int i,j,k,count;
        t_atom *itom,*jtom,*ktom;
@@ -833,9 +847,6 @@ printf("start pot force calc\n");
 
        /* get energy and force of every atom */
        for(i=0;i<count;i++) {
-printf("atom %d: %f\n",i,itom[i].r.x);
-printf("atom %d: %f\n",i,itom[i].v.x);
-printf("atom %d: %f\n",i,itom[i].f.x);
 
                /* reset force */
                v3_zero(&(itom[i].f));
@@ -844,95 +855,100 @@ printf("atom %d: %f\n",i,itom[i].f.x);
                if(itom[i].attr&ATOM_ATTR_1BP)
                        moldyn->func1b(moldyn,&(itom[i]));
 
+               if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
+                       continue;
+
                /* 2 body pair potential/force */
-               if(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)) {
        
-                       link_cell_neighbour_index(moldyn,
-                               (itom[i].r.x+moldyn->dim.x/2)/lc->x,
-                               (itom[i].r.y+moldyn->dim.y/2)/lc->y,
-                               (itom[i].r.z+moldyn->dim.z/2)/lc->z,
-                               neighbour_i);
+               link_cell_neighbour_index(moldyn,
+                                         (itom[i].r.x+moldyn->dim.x/2)/lc->x,
+                                         (itom[i].r.y+moldyn->dim.y/2)/lc->y,
+                                         (itom[i].r.z+moldyn->dim.z/2)/lc->z,
+                                         neighbour_i);
 
-                       dnlc=lc->dnlc;
+               dnlc=lc->dnlc;
 
-                       for(j=0;j<27;j++) {
+               for(j=0;j<27;j++) {
 
-                               this=&(neighbour_i[j]);
-                               list_reset(this);
+                       this=&(neighbour_i[j]);
+                       list_reset(this);
 
-                               if(this->start==NULL)
-                                       continue;
+                       if(this->start==NULL)
+                               continue;
 
-                               bc_ij=(j<dnlc)?0:1;
+                       bc_ij=(j<dnlc)?0:1;
 
-                               do {
-                                       jtom=this->current->data;
+                       do {
+                               jtom=this->current->data;
 
-                                       if(jtom==&(itom[i]))
-                                               continue;
+                               if(jtom==&(itom[i]))
+                                       continue;
 
-                                       if((jtom->attr&ATOM_ATTR_2BP)&
-                                          (itom[i].attr&ATOM_ATTR_2BP))
-                                               moldyn->func2b(moldyn,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
+                               if((jtom->attr&ATOM_ATTR_2BP)&
+                                  (itom[i].attr&ATOM_ATTR_2BP))
+                                       moldyn->func2b(moldyn,
+                                                      &(itom[i]),
+                                                      jtom,
+                                                      bc_ij);
 
-                                       /* 3 body potential/force */
+                               /* 3 body potential/force */
 
-                                       if(!(itom[i].attr&ATOM_ATTR_3BP)||
-                                          !(jtom->attr&ATOM_ATTR_3BP))
-                                               continue;
+                               if(!(itom[i].attr&ATOM_ATTR_3BP)||
+                                  !(jtom->attr&ATOM_ATTR_3BP))
+                                       continue;
 
-                                       /* copy the neighbour lists */
-                                       memcpy(neighbour_i2,neighbour_i,
-                                              27*sizeof(t_list));
+                               /* copy the neighbour lists */
+                               memcpy(neighbour_i2,neighbour_i,
+                                      27*sizeof(t_list));
 
-                                       /* get neighbours of i */
-                                       for(k=0;k<27;k++) {
+                               /* get neighbours of i */
+                               for(k=0;k<27;k++) {
 
-                                               that=&(neighbour_i2[k]);
-                                               list_reset(that);
+                                       that=&(neighbour_i2[k]);
+                                       list_reset(that);
                                        
-                                               if(that->start==NULL)
-                                                       continue;
+                                       if(that->start==NULL)
+                                               continue;
 
-                                               bc_ik=(k<dnlc)?0:1;
+                                       bc_ik=(k<dnlc)?0:1;
 
-                                               do {
+                                       do {
 
-                       ktom=that->current->data;
+                                               ktom=that->current->data;
 
-                       if(!(ktom->attr&ATOM_ATTR_3BP))
-                               continue;
+                                               if(!(ktom->attr&ATOM_ATTR_3BP))
+                                                       continue;
 
-                       if(ktom==jtom)
-                               continue;
+                                               if(ktom==jtom)
+                                                       continue;
 
-                       if(ktom==&(itom[i]))
-                               continue;
+                                               if(ktom==&(itom[i]))
+                                                       continue;
 
-                       moldyn->func3b(moldyn,&(itom[i]),jtom,ktom,bc_ik|bc_ij);
+                                               moldyn->func3b(moldyn,
+                                                              &(itom[i]),
+                                                              jtom,
+                                                              ktom,
+                                                              bc_ik|bc_ij);
 
-                                               } while(list_next(that)!=\
-                                                       L_NO_NEXT_ELEMENT);
+                                       } while(list_next(that)!=\
+                                               L_NO_NEXT_ELEMENT);
+
+                               }
 
-                                       }
-                                       
-                               } while(list_next(this)!=L_NO_NEXT_ELEMENT);
-               
                                /* 2bp post function */
                                if(moldyn->func2b_post) {
                                        moldyn->func2b_post(moldyn,
                                                            &(itom[i]),
                                                            jtom,bc_ij);
                                }
-
-                       }
+                                       
+                       } while(list_next(this)!=L_NO_NEXT_ELEMENT);
+               
                }
+
        }
 
-printf("end pot force calc\n");
        return 0;
 }
 
index a2223c0345472043779822c5aa80883c2665b464..72ff30bdfb00ab876cce6261360b827298c77043 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -92,6 +92,7 @@ typedef struct s_moldyn {
        double p;               /* actual pressure */
 
        /* pressure and temperature control (velocity/volume scaling) */
+       /* (in units of tau!) */
        unsigned char pt_scale; /* type of p and t scaling */
        double t_tc;            /* t berendsen control time constant */
        double p_tc;            /* p berendsen control time constant */
@@ -111,17 +112,15 @@ typedef struct s_moldyn {
        double energy;          /* potential energy */
        double ekin;            /* kinetic energy */
 
-       t_visual vis;           /* visualization/log/save interface structure */
-       u8 lvstat;              /* log & vis properties */
+       char vlsdir[128];       /* visualization/log/save directory */
+       t_visual vis;           /* visualization interface structure */
+       u8 vlsprop;             /* log/vis/save properties */
        unsigned int ewrite;    /* how often to log energy */
        int efd;                /* fd for energy log */
        unsigned int mwrite;    /* how often to log momentum */
        int mfd;                /* fd for momentum log */
        unsigned int vwrite;    /* how often to visualize atom information */
-       char vfb[64];           /* visualization file name base */
-       //void *visual;         /* pointer (hack!) */
        unsigned int swrite;    /* how often to create a save file */
-       char sfb[64];           /* visualization file name base */
 
        u8 status;              /* general moldyn properties */
 
@@ -359,7 +358,8 @@ int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params);
 int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params);
 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params);
 
-int moldyn_set_log(t_moldyn *moldyn,u8 type,char *fb,int timer);
+int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
+int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
 int moldyn_log_shutdown(t_moldyn *moldyn);
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
diff --git a/ppm2avi b/ppm2avi
index be4c6dd9555ab85422e912e9c3c0a9890d0e0e9d..4a630a73fdb8b11fbe448374fa16bd88a9c79dc0 100755 (executable)
--- a/ppm2avi
+++ b/ppm2avi
@@ -1,5 +1,5 @@
-for i in saves/*.ppm; do
-       convert $i saves/$(basename $i ppm)png
+for i in $1/*.ppm; do
+       convert $i $1/$(basename $i ppm)png
 done
 
-mencoder mf://saves/*.png -ovc copy -o video/md.avi
+mencoder mf://$1/*.png -ovc copy -o video/md.avi
diff --git a/run b/run
index 264a4aca518259049bfdeaea37f1c26e3d4885f4..af5d425b9a19931c2f44bb0d011ee4269859e018 100755 (executable)
--- a/run
+++ b/run
@@ -1,10 +1,10 @@
 mkdir -p saves video
-./clean
+./clean $1
 ./sic
 if [ "$?" == "0" ]; then
-       ./perms
+       #./perms
        if [ "$1" ] ; then
-               rasmol -nodisplay < $1.scr > /dev/null 2>&1
-               ./ppm2avi
+               rasmol -nodisplay < $1/visualize.scr > /dev/null 2>&1
+               ./ppm2avi $1
        fi
 fi
diff --git a/sic.c b/sic.c
index 8c0e648cf544ad454d7a45d9421314213cc3a31a..1e29a85d9647341d316def20e9f1a097046d41f3 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -113,13 +113,13 @@ int main(int argc,char **argv) {
        //               0,5,5,5);
 
        /* testing configuration */
-       r.x=-0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
-       r.y=0;                  v.y=0;
-       r.z=0;                  v.z=0;
+       r.x=2.7/2;      v.x=0;
+       r.y=0;          v.y=0;
+       r.z=0;          v.z=0;
        add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
-       r.x=+0.55*0.25*sqrt(3.0)*LC_SI; v.x=0;
-       r.y=0;                  v.y=0;
-       r.z=0;                  v.z=0;
+       r.x=-2.7/2;     v.x=0;
+       r.y=0;          v.y=0;
+       r.z=0;          v.z=0;
        add_atom(&md,SI,M_SI,0,ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,&r,&v);
 
        /* setting a nearest neighbour distance for the moldyn checks */
@@ -132,7 +132,7 @@ int main(int argc,char **argv) {
 
        /* set p/t scaling */
        printf("[sic] set p/t scaling\n");
-       set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100*tau);
+       set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
        
        /* initial thermal fluctuations of particles (in equilibrium) */
        printf("[sic] thermal init\n");
@@ -140,12 +140,13 @@ int main(int argc,char **argv) {
 
        /* create the simulation schedule */
        printf("[sic] adding schedule\n");
-       moldyn_add_schedule(&md,100,1.0);
+       moldyn_add_schedule(&md,1000,1.0);
 
        /* activate logging */
        printf("[sic] activate logging\n");
-       moldyn_set_log(&md,LOG_TOTAL_ENERGY,"saves/test-energy",1);
-       moldyn_set_log(&md,VISUAL_STEP,"saves/test-visual",1);
+       moldyn_set_log_dir(&md,"saves/test");
+       moldyn_set_log(&md,LOG_TOTAL_ENERGY,10);
+       moldyn_set_log(&md,VISUAL_STEP,10);
 
        /*
         * let's do the actual md algorithm now
index d03903bc9283c8fdefd818b51fe7bc6a2078474f..a6fda4f747082c714c085b299aaf62f24466031b 100644 (file)
@@ -46,9 +46,9 @@ int visual_init(t_visual *v,char *filebase) {
 
        strncpy(v->fb,filebase,128);
        memset(file,0,128+8);
-       sprintf(file,"%s.scr",v->fb);
+       sprintf(file,"%s/visualize.scr",v->fb);
 
-       v->fd=open(file,O_WRONLY|O_CREAT|O_TRUNC);
+       v->fd=open(file,O_WRONLY|O_CREAT|O_EXCL,S_IRUSR|S_IWUSR);
        if(v->fd<0) {
                perror("open visual fd");
                return -1;
@@ -78,8 +78,8 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dim.y=v->dim.y;
        dim.z=v->dim.z;
 
-       sprintf(file,"%s-%07f.xyz",v->fb,time);
-       fd=open(file,O_WRONLY|O_CREAT|O_TRUNC);
+       sprintf(file,"%s/visualize_%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;
@@ -94,13 +94,13 @@ int visual_atoms(t_visual *v,double time,t_atom *atom,int n) {
        dprintf(v->fd,"set specular on\n");
        dprintf(v->fd,"set boundbox on\n");
        //dprintf(v->fd,"label\n");
-       sprintf(file,"%s-%.15f.ppm",v->fb,time);
+       sprintf(file,"%s/visualize_%07.f.ppm",v->fb,time);
        dprintf(v->fd,"write ppm %s\n",file);
        dprintf(v->fd,"zap\n");
 
        /* write the actual data file */
        dprintf(fd,"%d\n",(dim.x==0)?n:n+8);
-       dprintf(fd,"atoms at time %.15f\n",time);
+       dprintf(fd,"atoms at time %07.f\n",time);
        for(i=0;i<n;i++)
                dprintf(fd,"%s %f %f %f\n",pse[atom[i].element],
                                           atom[i].r.x,