safety checkin
[physik/posic.git] / moldyn.c
index 4ddc42e..7ea5d50 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,13 +404,14 @@ 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;
+                       scale=1.0+(scale-1.0)/moldyn->t_tc;
        scale=sqrt(scale);
 
        /* velocity scaling */
-       for(i=0;i<moldyn->count;i++)
+       for(i=0;i<moldyn->count;i++) {
                if((equi_init&TRUE)||(atom[i].attr&ATOM_ATTR_HB))
                        v3_scale(&(atom[i].v),&(atom[i].v),scale);
+       }
 
        return 0;
 }
@@ -447,12 +468,6 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist) {
 
        /* nn_dist is the nearest neighbour distance */
 
-       if(moldyn->t==5.0) {
-               printf("[moldyn] i do not estimate timesteps below %f K!\n",
-                      MOLDYN_CRITICAL_EST_TEMP);
-               return 23.42;
-       }
-
        tau=(0.05*nn_dist*moldyn->atom[0].mass)/sqrt(3.0*K_BOLTZMANN*moldyn->t);
 
        return tau;     
@@ -645,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);
@@ -703,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));
@@ -719,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));
@@ -747,6 +758,9 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        }
                }
 
+               /* increase absolute time */
+               moldyn->time+=moldyn->tau;
+
        }
 
                /* check for hooks */
@@ -781,15 +795,12 @@ int velocity_verlet(t_moldyn *moldyn) {
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
                v3_scale(&delta,&(atom[i].f),0.5*tau_square/atom[i].mass);
                v3_add(&(atom[i].r),&(atom[i].r),&delta);
-//if(i==0) printf("und dann %.20f\n",atom[i].r.x*1e10);
                check_per_bound(moldyn,&(atom[i].r));
-//if(i==0) printf("und danach %.20f\n",atom[i].r.x);
 
                /* velocities */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
        }
-moldyn_bc_check(moldyn);
 
        /* neighbour list update */
        link_cell_update(moldyn);
@@ -844,93 +855,98 @@ int potential_force_calc(t_moldyn *moldyn) {
                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("DEBUG: %.15f \n",itom[i].f.x);
+
        }
 
        return 0;
@@ -1053,8 +1069,8 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
        p->mu_m=0.5*(p->mu[0]+p->mu[1]);
 
        printf("[moldyn] tersoff mult parameter info:\n");
-       printf("  S (m)  | %.12f | %.12f | %.12f\n",p->S[0],p->S[1],p->Smixed);
-       printf("  R (m)  | %.12f | %.12f | %.12f\n",p->R[0],p->R[1],p->Rmixed);
+       printf("  S (A)  | %f | %f | %f\n",p->S[0],p->S[1],p->Smixed);
+       printf("  R (A)  | %f | %f | %f\n",p->R[0],p->R[1],p->Rmixed);
        printf("  A (eV) | %f | %f | %f\n",p->A[0]/EV,p->A[1]/EV,p->Amixed/EV);
        printf("  B (eV) | %f | %f | %f\n",p->B[0]/EV,p->B[1]/EV,p->Bmixed/EV);
        printf("  lambda | %f | %f | %f\n",p->lambda[0],p->lambda[1],