extended debugging, implemented fixed position attribute

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index ec96712..80cf572 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -595,7 +595,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        int count;
        
        atom=moldyn->atom;
-       count=(moldyn->count)++;
+       count=(moldyn->count)++;        // asshole style!
 
        ptr=realloc(atom,(count+1)*sizeof(t_atom));
        if(!ptr) {
@@ -605,6 +605,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        moldyn->atom=ptr;
 
        atom=moldyn->atom;
+
+       /* initialize new atom */
+       memset(&(atom[count]),0,sizeof(t_atom));
        atom[count].r=*r;
        atom[count].v=*v;
        atom[count].element=element;
@@ -956,6 +959,8 @@ double pressure_calc(t_moldyn *moldyn) {
 
 int average_reset(t_moldyn *moldyn) {
 
+       printf("[moldyn] average reset\n");
+
        /* update skip value */
        moldyn->avg_skip=moldyn->total_steps;
 
@@ -1369,7 +1374,7 @@ int link_cell_update(t_moldyn *moldyn) {
                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)
@@ -1419,6 +1424,10 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,
        count2=27;
        a=nx*ny;
 
+       if(i>=nx||j>=ny||k>=nz)
+               printf("[moldyn] WARNING: lcni %d/%d %d/%d %d/%d\n",
+                      i,nx,j,ny,k,nz);
+
        cell[0]=lc->subcell[i+j*nx+k*a];
        for(ci=-1;ci<=1;ci++) {
                bx=0;
@@ -1569,14 +1578,14 @@ int moldyn_integrate(t_moldyn *moldyn) {
 
        /* some stupid checks before we actually start calculating bullshit */
        if(moldyn->cutoff>0.5*moldyn->dim.x)
-               printf("[moldyn] warning: cutoff > 0.5 x dim.x\n");
+               printf("[moldyn] WARNING: cutoff > 0.5 x dim.x\n");
        if(moldyn->cutoff>0.5*moldyn->dim.y)
-               printf("[moldyn] warning: cutoff > 0.5 x dim.y\n");
+               printf("[moldyn] WARNING: cutoff > 0.5 x dim.y\n");
        if(moldyn->cutoff>0.5*moldyn->dim.z)
-               printf("[moldyn] warning: cutoff > 0.5 x dim.z\n");
+               printf("[moldyn] WARNING: cutoff > 0.5 x dim.z\n");
        ds=0.5*atom[0].f.x*moldyn->tau_square/atom[0].mass;
        if(ds>0.05*moldyn->nnd)
-               printf("[moldyn] warning: forces too high / tau too small!\n");
+               printf("[moldyn] WARNING: forces too high / tau too small!\n");
 
        /* zero absolute time */
        moldyn->time=0.0;
@@ -1729,6 +1738,9 @@ int velocity_verlet(t_moldyn *moldyn) {
        tau_square=moldyn->tau_square;
 
        for(i=0;i<count;i++) {
+               /* check whether fixed atom */
+               if(atom[i].attr&ATOM_ATTR_FP)
+                       continue;
                /* new positions */
                h=0.5/atom[i].mass;
                v3_scale(&delta,&(atom[i].v),tau);
@@ -1752,6 +1764,9 @@ int velocity_verlet(t_moldyn *moldyn) {
        potential_force_calc(moldyn);
 
        for(i=0;i<count;i++) {
+               /* check whether fixed atom */
+               if(atom[i].attr&ATOM_ATTR_FP)
+                       continue;
                /* again velocities [actually v(t+tau)] */
                v3_scale(&delta,&(atom[i].f),0.5*tau/atom[i].mass);
                v3_add(&(atom[i].v),&(atom[i].v),&delta);
@@ -2079,20 +2094,27 @@ int potential_force_calc(t_moldyn *moldyn) {
        //printf("\nATOM 0: %f %f %f\n\n",itom->f.x,itom->f.y,itom->f.z);
        if(moldyn->time>DSTART&&moldyn->time<DEND) {
                printf("force:\n");
-               printf("  x: %0.40f\n",moldyn->atom[5832].f.x);
-               printf("  y: %0.40f\n",moldyn->atom[5832].f.y);
-               printf("  z: %0.40f\n",moldyn->atom[5832].f.z);
+               printf("  x: %0.40f\n",moldyn->atom[DATOM].f.x);
+               printf("  y: %0.40f\n",moldyn->atom[DATOM].f.y);
+               printf("  z: %0.40f\n",moldyn->atom[DATOM].f.z);
        }
 #endif
 
-       /* calculate global virial */
+       /* some postprocessing */
        for(i=0;i<count;i++) {
-               moldyn->gvir.xx+=moldyn->atom[i].r.x*moldyn->atom[i].f.x;
-               moldyn->gvir.yy+=moldyn->atom[i].r.y*moldyn->atom[i].f.y;
-               moldyn->gvir.zz+=moldyn->atom[i].r.z*moldyn->atom[i].f.z;
-               moldyn->gvir.xy+=moldyn->atom[i].r.y*moldyn->atom[i].f.x;
-               moldyn->gvir.xz+=moldyn->atom[i].r.z*moldyn->atom[i].f.x;
-               moldyn->gvir.yz+=moldyn->atom[i].r.z*moldyn->atom[i].f.y;
+               /* calculate global virial */
+               moldyn->gvir.xx+=itom[i].r.x*itom[i].f.x;
+               moldyn->gvir.yy+=itom[i].r.y*itom[i].f.y;
+               moldyn->gvir.zz+=itom[i].r.z*itom[i].f.z;
+               moldyn->gvir.xy+=itom[i].r.y*itom[i].f.x;
+               moldyn->gvir.xz+=itom[i].r.z*itom[i].f.x;
+               moldyn->gvir.yz+=itom[i].r.z*itom[i].f.y;
+
+               /* check forces regarding the given timestep */
+               if(v3_norm(&(itom[i].f))>\
+                  0.1*moldyn->nnd*itom[i].mass/moldyn->tau_square)
+                       printf("[moldyn] WARNING: pfc (high force: atom %d)\n",
+                              i);
        }
 
        return 0;
@@ -2207,6 +2229,8 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 
        int fd;
        int cnt,size;
+       int fsize;
+       int corr;
 
        fd=open(file,O_RDONLY);
        if(fd<0) {
@@ -2214,6 +2238,9 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
                return fd;
        }
 
+       fsize=lseek(fd,0,SEEK_END);
+       lseek(fd,0,SEEK_SET);
+
        size=sizeof(t_moldyn);
 
        while(size) {
@@ -2227,6 +2254,19 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
 
        size=moldyn->count*sizeof(t_atom);
 
+       /* correcting possible atom data offset */
+       corr=0;
+       if(fsize!=sizeof(t_moldyn)+size) {
+               corr=fsize-sizeof(t_moldyn)-size;
+               printf("[moldyn] WARNING: lsf (illegal file size)\n");
+               printf("  moifying offset:\n");
+               printf("  - current pos: %d\n",sizeof(t_moldyn));
+               printf("  - atom size: %d\n",size);
+               printf("  - file size: %d\n",fsize);
+               printf("  => correction: %d\n",corr);
+               lseek(fd,corr,SEEK_CUR);
+       }
+
        moldyn->atom=(t_atom *)malloc(size);
        if(moldyn->atom==NULL) {
                perror("[moldyn] load save file malloc (atoms)");
@@ -2308,6 +2348,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        slots=moldyn->cutoff/dr;
        o=2*slots;
 
+       if(slots*dr<=moldyn->cutoff)
+               printf("[moldyn] WARNING: pcc (low #slots)\n");
+
        printf("[moldyn] pair correlation calc info:\n");
        printf("  time: %f\n",moldyn->time);
        printf("  count: %d\n",moldyn->count);
@@ -2389,8 +2432,9 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
                                /* should never happen but it does 8) -
                                 * related to -ffloat-store problem! */
                                if(s>=slots) {
-                                       printf("[moldyn] WARNING pcc (%d/%d)\n",
+                                       printf("[moldyn] WARNING: pcc (%d/%d)",
                                               s,slots);
+                                       printf("\n");
                                        s=slots-1;
                                }