moldyn->vis.dim.z=z;
}
+ printf("[moldyn] dimensions in A:\n");
+ printf(" x: %f\n",moldyn->dim.x);
+ printf(" y: %f\n",moldyn->dim.y);
+ printf(" z: %f\n",moldyn->dim.z);
+ printf(" visualize simulation box: %s\n",visualize?"on":"off");
+
return 0;
}
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;
}
count=a*b*c;
+ /* how many atoms do we expect */
if(type==FCC) count*=4;
-
if(type==DIAMOND) count*=8;
+ /* allocate space for atoms */
moldyn->atom=malloc(count*sizeof(t_atom));
if(moldyn->atom==NULL) {
perror("malloc (atoms)");
/* debug */
if(ret!=count) {
- printf("ok, there is something wrong ...\n");
- printf("calculated -> %d atoms\n",count);
- printf("created -> %d atoms\n",ret);
+ printf("[moldyn] creating lattice failed\n");
+ printf(" amount of atoms\n");
+ printf(" - expected: %d\n",count);
+ printf(" - created: %d\n",ret);
return -1;
}
check_per_bound(moldyn,&(moldyn->atom[count].r));
}
-
return ret;
}
/* (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 */
v3_add(&(atom[i].v),&(atom[i].v),&delta);
}
+moldyn_bc_check(moldyn);
/* neighbour list update */
link_cell_update(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;
/* 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;
}
dim=&(moldyn->dim);
- x=0.5*dim->x;
- y=0.5*dim->y;
- z=0.5*dim->z;
+ x=dim->x/2;
+ y=dim->y/2;
+ z=dim->z/2;
if(moldyn->status&MOLDYN_STAT_PBX) {
if(a->x>=x) a->x-=dim->x;
if(a->z>=z) a->z-=dim->z;
else if(-a->z>z) a->z+=dim->z;
}
+printf("%f %f %f\n",a->x,x,a->x/x);
return 0;
}
/* f_a, df_a calc (again, same for ij and ji) | save for later use! */
exchange->f_a=-B*exp(-mu*d_ij);
- exchange->df_a=-mu*exchange->f_a/d_ij;
+ exchange->df_a=mu*exchange->f_a/d_ij;
/* f_c, df_c calc (again, same for ij and ji) */
if(d_ij<R) {
v3_scale(&force,&force,-0.5);
/* add force */
- v3_sub(&(ai->f),&(ai->f),&force);
+ v3_add(&(ai->f),&(ai->f),&force);
return 0;
}
for(i=0;i<moldyn->count;i++) {
if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2)
printf("FATAL: atom %d: x: %.20f (%.20f)\n",
- i,atom[i].r.x*1e10,dim->x/2*1e10);
+ i,atom[i].r.x,dim->x/2);
if(atom[i].r.y>=dim->y/2||-atom[i].r.y>dim->y/2)
printf("FATAL: atom %d: y: %.20f (%.20f)\n",
- i,atom[i].r.y*1e10,dim->y/2*1e10);
+ i,atom[i].r.y,dim->y/2);
if(atom[i].r.z>=dim->z/2||-atom[i].r.z>dim->z/2)
printf("FATAL: atom %d: z: %.20f (%.20f)\n",
- i,atom[i].r.z*1e10,dim->z/2*1e10);
+ i,atom[i].r.z,dim->z/2);
}
return 0;
projects / physik / posic.git / blobdiff