int average_and_fluctuation_calc(t_moldyn *moldyn) {
+ int denom;
+
if(moldyn->total_steps<moldyn->avg_skip)
return 0;
- int denom=moldyn->total_steps+1-moldyn->avg_skip;
+ denom=moldyn->total_steps+1-moldyn->avg_skip;
/* assume up to date energies, temperature, pressure etc */
double thermodynamic_pressure_calc(t_moldyn *moldyn) {
- t_3dvec dim,*tp;
+ t_3dvec dim;
+ //t_3dvec *tp;
double u_up,u_down,dv;
double scale,p;
t_atom *store;
/* linked list / cell method */
-#ifdef STATIC_LISTS
-
int link_cell_init(t_moldyn *moldyn,u8 vol) {
t_linkcell *lc;
int i;
- int *foo;
lc=&(moldyn->lc);
lc->y=moldyn->dim.y/lc->ny;
lc->nz=moldyn->dim.z/moldyn->cutoff;
lc->z=moldyn->dim.z/lc->nz;
-
lc->cells=lc->nx*lc->ny*lc->nz;
+
+#ifdef STATIC_LISTS
lc->subcell=malloc(lc->cells*sizeof(int*));
+#else
+ lc->subcell=malloc(lc->cells*sizeof(t_list));
+#endif
+
+ if(lc->subcell==NULL) {
+ perror("[moldyn] cell init (malloc)");
+ return -1;
+ }
if(lc->cells<27)
printf("[moldyn] FATAL: less then 27 subcells!\n");
if(vol) {
+#ifdef STATIC_LISTS
printf("[moldyn] initializing 'static' linked cells (%d)\n",
lc->cells);
+#else
+ printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
+ lc->cells);
+#endif
printf(" x: %d x %f A\n",lc->nx,lc->x);
printf(" y: %d x %f A\n",lc->ny,lc->y);
printf(" z: %d x %f A\n",lc->nz,lc->z);
}
+#ifdef STATIC_LISTS
/* list init */
for(i=0;i<lc->cells;i++) {
lc->subcell[i]=malloc((MAX_ATOMS_PER_LIST+1)*sizeof(int));
perror("[moldyn] list init (malloc)");
return -1;
}
-//if(i==0) printf(" --- add one here! %d %p %p ----\n",i,lc->subcell,lc->subcell[0]);
+ /*
+ if(i==0)
+ printf(" ---> %d malloc %p (%p)\n",
+ i,lc->subcell[0],lc->subcell);
+ */
}
+#else
+ for(i=0;i<lc->cells;i++)
+ list_init_f(&(lc->subcell[i]));
+#endif
/* update the list */
link_cell_update(moldyn);
int nx,ny;
t_atom *atom;
t_linkcell *lc;
+#ifdef STATIC_LISTS
int p;
+#endif
atom=moldyn->atom;
lc=&(moldyn->lc);
ny=lc->ny;
for(i=0;i<lc->cells;i++)
+#ifdef STATIC_LISTS
memset(lc->subcell[i],0,(MAX_ATOMS_PER_LIST+1)*sizeof(int));
+#else
+ list_destroy_f(&(lc->subcell[i]));
+#endif
for(count=0;count<moldyn->count;count++) {
i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
+#ifdef STATIC_LISTS
p=0;
while(lc->subcell[i+j*nx+k*nx*ny][p]!=0)
p++;
}
lc->subcell[i+j*nx+k*nx*ny][p]=count;
+#else
+ list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
+ &(atom[count]));
+ /*
+ if(j==0&&k==0)
+ printf(" ---> %d %d malloc %p (%p)\n",
+ i,count,lc->subcell[i].current,lc->subcell);
+ */
+#endif
}
return 0;
}
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell) {
+int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
+#ifdef STATIC_LISTS
+ int **cell
+#else
+ t_list *cell
+#endif
+ ) {
t_linkcell *lc;
int a;
int i;
t_linkcell *lc;
- int *foo;
lc=&(moldyn->lc);
- for(i=0;i<lc->cells;i++)
-{
-//printf(" --- free %p , %d\n",lc->subcell[i],i);
+ for(i=0;i<lc->cells;i++) {
+#ifdef STATIC_LISTS
free(lc->subcell[i]);
-}
-
- free(lc->subcell);
-
- return 0;
-}
-
#else
-
-int link_cell_init(t_moldyn *moldyn,u8 vol) {
-
- t_linkcell *lc;
- int i;
-
- lc=&(moldyn->lc);
-
- /* partitioning the md cell */
- lc->nx=moldyn->dim.x/moldyn->cutoff;
- lc->x=moldyn->dim.x/lc->nx;
- lc->ny=moldyn->dim.y/moldyn->cutoff;
- lc->y=moldyn->dim.y/lc->ny;
- lc->nz=moldyn->dim.z/moldyn->cutoff;
- lc->z=moldyn->dim.z/lc->nz;
-
- lc->cells=lc->nx*lc->ny*lc->nz;
- lc->subcell=malloc(lc->cells*sizeof(t_list));
-
- if(lc->cells<27)
- printf("[moldyn] FATAL: less then 27 subcells!\n");
-
- if(vol) {
- printf("[moldyn] initializing 'dynamic' linked cells (%d)\n",
- lc->cells);
- printf(" x: %d x %f A\n",lc->nx,lc->x);
- printf(" y: %d x %f A\n",lc->ny,lc->y);
- printf(" z: %d x %f A\n",lc->nz,lc->z);
- }
-
- for(i=0;i<lc->cells;i++)
- list_init_f(&(lc->subcell[i]));
-
- link_cell_update(moldyn);
-
- return 0;
-}
-
-int link_cell_update(t_moldyn *moldyn) {
-
- int count,i,j,k;
- int nx,ny;
- t_atom *atom;
- t_linkcell *lc;
- double x,y,z;
-
- atom=moldyn->atom;
- lc=&(moldyn->lc);
-
- nx=lc->nx;
- ny=lc->ny;
-
- x=moldyn->dim.x/2;
- y=moldyn->dim.y/2;
- z=moldyn->dim.z/2;
-
- for(i=0;i<lc->cells;i++)
+ //printf(" ---> %d free %p\n",i,lc->subcell[i].start);
list_destroy_f(&(lc->subcell[i]));
-
- for(count=0;count<moldyn->count;count++) {
- i=((atom[count].r.x+(moldyn->dim.x/2))/lc->x);
- j=((atom[count].r.y+(moldyn->dim.y/2))/lc->y);
- k=((atom[count].r.z+(moldyn->dim.z/2))/lc->z);
- list_add_immediate_f(&(lc->subcell[i+j*nx+k*nx*ny]),
- &(atom[count]));
-//if(i==0&&j==0&&k==0) printf(" --- add one here! %d %p ----\n",count,lc->subcell[0].current);
- }
-
- return 0;
-}
-
-int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) {
-
- t_linkcell *lc;
- int a;
- int count1,count2;
- int ci,cj,ck;
- int nx,ny,nz;
- int x,y,z;
- u8 bx,by,bz;
-
- lc=&(moldyn->lc);
- nx=lc->nx;
- ny=lc->ny;
- nz=lc->nz;
- count1=1;
- count2=27;
- a=nx*ny;
-
- cell[0]=lc->subcell[i+j*nx+k*a];
- for(ci=-1;ci<=1;ci++) {
- bx=0;
- x=i+ci;
- if((x<0)||(x>=nx)) {
- x=(x+nx)%nx;
- bx=1;
- }
- for(cj=-1;cj<=1;cj++) {
- by=0;
- y=j+cj;
- if((y<0)||(y>=ny)) {
- y=(y+ny)%ny;
- by=1;
- }
- for(ck=-1;ck<=1;ck++) {
- bz=0;
- z=k+ck;
- if((z<0)||(z>=nz)) {
- z=(z+nz)%nz;
- bz=1;
- }
- if(!(ci|cj|ck)) continue;
- if(bx|by|bz) {
- cell[--count2]=lc->subcell[x+y*nx+z*a];
- }
- else {
- cell[count1++]=lc->subcell[x+y*nx+z*a];
- }
- }
- }
+#endif
}
- lc->dnlc=count1;
-
- return count1;
-}
-
-int link_cell_shutdown(t_moldyn *moldyn) {
-
- int i;
- t_linkcell *lc;
-
- lc=&(moldyn->lc);
-
-printf("FOO:\n");
- for(i=0;i<lc->nx*lc->ny*lc->nz;i++) {
-printf(" %d\n",i);
- list_destroy_f(&(moldyn->lc.subcell[i]));
-printf(" %d!\n",i);
-}
-
free(lc->subcell);
return 0;
}
-#endif
-
int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau) {
int count;
/* get current time */
gettimeofday(&t2,NULL);
- printf("\rsched:%d, steps:%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
- sched->count,i,
- moldyn->t,moldyn->t_avg,
- moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
- moldyn->volume,
- (int)(t2.tv_sec-t1.tv_sec));
- fflush(stdout);
+printf("\rsched:%d, steps:%d/%d, T:%3.1f/%3.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+ sched->count,i,moldyn->total_steps,
+ moldyn->t,moldyn->t_avg,
+ moldyn->p_avg/BAR,moldyn->gp_avg/BAR,
+ moldyn->volume,
+ (int)(t2.tv_sec-t1.tv_sec));
+
+ fflush(stdout);
/* copy over time */
t1=t2;
/* check for hooks */
if(sched->hook) {
- printf("\n ## schedule hook %d/%d start ##\n",
- sched->count+1,sched->total_sched-1);
+ printf("\n ## schedule hook %d start ##\n",
+ sched->count);
sched->hook(moldyn,sched->hook_params);
printf(" ## schedule hook end ##\n");
}
}
size=sizeof(t_moldyn);
- cnt=read(fd,moldyn,size);
- if(cnt!=size) {
- perror("[moldyn] load save file read (moldyn)");
- return cnt;
+
+ while(size) {
+ cnt=read(fd,moldyn,size);
+ if(cnt<0) {
+ perror("[moldyn] load save file read (moldyn)");
+ return cnt;
+ }
+ size-=cnt;
}
size=moldyn->count*sizeof(t_atom);
return -1;
}
- cnt=read(fd,moldyn->atom,size);
- if(cnt!=size) {
- perror("[moldyn] load save file read (atoms)");
- return cnt;
+ while(size) {
+ cnt=read(fd,moldyn->atom,size);
+ if(cnt<0) {
+ perror("[moldyn] load save file read (atoms)");
+ return cnt;
+ }
+ size-=cnt;
}
- // hooks
+ // hooks etc ...
return 0;
}
t_list *this;
unsigned char bc;
t_3dvec dist;
- double d,norm;
+ double d;
+ //double norm;
int o,s;
unsigned char ibrand;
lc=&(moldyn->lc);
- slots=(int)(moldyn->cutoff/dr);
+ slots=moldyn->cutoff/dr;
o=2*slots;
printf("[moldyn] pair correlation calc info:\n");
jtom=this->current->data;
#endif
-
- if(jtom==&(itom[i]))
- continue;
-
- /* only count pairs once */
- if(itom[i].tag>jtom->tag)
+ /* only count pairs once,
+ * skip same atoms */
+ if(itom[i].tag>=jtom->tag)
continue;
/*
if(bc) check_per_bound(moldyn,&dist);
d=v3_absolute_square(&dist);
- /* ignore if greater cutoff */
- if(d>moldyn->cutoff_square)
+ /* ignore if greater or equal cutoff */
+ if(d>=moldyn->cutoff_square)
continue;
/* fill the slots */
d=sqrt(d);
s=(int)(d/dr);
+ /* should never happen but it does 8) -
+ * related to -ffloat-store problem! */
+ if(s>=slots) s=slots-1;
+
if(ibrand!=jtom->brand) {
/* mixed */
stat[s]+=1;
/* type b - type b bonds */
stat[s+o]+=1;
}
-
#ifdef STATIC_LISTS
}
#else
projects / physik / posic.git / blobdiff