From 1d83ceb2ce2ff5150fd079f1066b7f583e38c8f4 Mon Sep 17 00:00:00 2001 From: hackbard Date: Tue, 12 Dec 2006 11:50:04 +0000 Subject: [PATCH] safety checkin --- clean | 2 +- moldyn.c | 184 ++++++++++++++++++++++++++---------------------- moldyn.h | 12 ++-- ppm2avi | 6 +- run | 8 +-- sic.c | 21 +++--- visual/visual.c | 12 ++-- 7 files changed, 131 insertions(+), 114 deletions(-) diff --git a/clean b/clean index ff74c9b..ca57788 100755 --- a/clean +++ b/clean @@ -1 +1 @@ -rm -f saves/* video/*.ppm video/*.png > /dev/null 2>&1 +rm -f $1/* > /dev/null 2>&1 diff --git a/moldyn.c b/moldyn.c index 9d9901c..7ea5d50 100644 --- 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;icount;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;ifunc1b(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=(jcurrent->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=(kcurrent->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; } diff --git a/moldyn.h b/moldyn.h index a2223c0..72ff30b 100644 --- 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 be4c6dd..4a630a7 100755 --- 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 264a4ac..af5d425 100755 --- 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 8c0e648..1e29a85 100644 --- 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 diff --git a/visual/visual.c b/visual/visual.c index d03903b..a6fda4f 100644 --- a/visual/visual.c +++ b/visual/visual.c @@ -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