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;
}
/* (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 */
}
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;
t_moldyn_schedule *schedule;
t_atom *atom;
int fd;
- char fb[128];
+ char dir[128];
double ds;
schedule=&(moldyn->schedule);
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));
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));
}
}
+ /* increase absolute time */
+ moldyn->time+=moldyn->tau;
+
}
/* check for hooks */
/* 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;
/* 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));
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;
}
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 */
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 */
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,
// 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 */
/* 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");
/* 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
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;
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;
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,
projects / physik / posic.git / commitdiff