regular checkin of logfile
[physik/posic.git] / moldyn.c
index d7d89bc..bc4fdff 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -23,7 +23,7 @@
 #include <omp.h>
 #endif
 
-#ifdef PTHREADS
+#if defined PTHREADS || defined VISUAL_THREAD
 #include <pthread.h>
 #endif
 
 #endif
 
 /* pse */
+#define PSE_MASS
 #define PSE_NAME
 #define PSE_COL
 #include "pse.h"
+#undef PSE_MASS
 #undef PSE_NAME
 #undef PSE_COL
 
+#ifdef PTHREADS
+/* global mutexes */
+pthread_mutex_t *amutex;
+pthread_mutex_t emutex;
+#endif
+
 /*
  * the moldyn functions
  */
@@ -66,13 +74,28 @@ int moldyn_init(t_moldyn *moldyn,int argc,char **argv) {
        rand_init(&(moldyn->random),NULL,1);
        moldyn->random.status|=RAND_STAT_VERBOSE;
 
+#ifdef PTHREADS
+       pthread_mutex_init(&emutex,NULL);
+#endif
+
        return 0;
 }
 
 int moldyn_shutdown(t_moldyn *moldyn) {
 
+#ifdef PTHREADS
+       int i;
+#endif
+
        printf("[moldyn] shutdown\n");
 
+#ifdef PTHREADS
+       for(i=0;i<moldyn->count;i++)
+               pthread_mutex_destroy(&(amutex[i]));
+       free(amutex);
+       pthread_mutex_destroy(&emutex);
+#endif
+
        moldyn_log_shutdown(moldyn);
        link_cell_shutdown(moldyn);
        rand_close(&(moldyn->random));
@@ -491,8 +514,9 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  * creating lattice functions
  */
 
-int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin) {
+int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,
+                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
+                   t_part_params *p_params,t_defect_params *d_params) {
 
        int new,count;
        int ret;
@@ -500,6 +524,9 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        void *ptr;
        t_atom *atom;
        char name[16];
+#ifdef PTHREADS
+       pthread_mutex_t *mutex;
+#endif
 
        new=a*b*c;
        count=moldyn->count;
@@ -513,6 +540,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        if(type==FCC) new*=4;
        if(type==DIAMOND) new*=8;
 
+       /* defects */
+       if(d_params->type) {
+               switch(d_params->stype) {
+                       case DEFECT_STYPE_DB_X:
+                       case DEFECT_STYPE_DB_Y:
+                       case DEFECT_STYPE_DB_Z:
+                       case DEFECT_STYPE_DB_R:
+                               new*=2;
+                               break;
+                       default:
+                               printf("[moldyn] WARNING: cl unknown defect\n");
+                               break;
+               }
+       }
+
        /* allocate space for atoms */
        ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
        if(!ptr) {
@@ -522,6 +564,16 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        moldyn->atom=ptr;
        atom=&(moldyn->atom[count]);
 
+#ifdef PTHREADS
+       ptr=realloc(amutex,(count+new)*sizeof(pthread_mutex_t));
+       if(!ptr) {
+               perror("[moldyn] mutex realloc (add atom)");
+               return -1;
+       }
+       amutex=ptr;
+       mutex=&(amutex[count]);
+#endif
+
        /* no atoms on the boundaries (only reason: it looks better!) */
        if(!origin) {
                orig.x=0.5*lc;
@@ -537,21 +589,21 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        switch(type) {
                case CUBIC:
                        set_nn_dist(moldyn,lc);
-                       ret=cubic_init(a,b,c,lc,atom,&orig);
+                       ret=cubic_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"cubic");
                        break;
                case FCC:
                        if(!origin)
                                v3_scale(&orig,&orig,0.5);
                        set_nn_dist(moldyn,0.5*sqrt(2.0)*lc);
-                       ret=fcc_init(a,b,c,lc,atom,&orig);
+                       ret=fcc_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"fcc");
                        break;
                case DIAMOND:
                        if(!origin)
                                v3_scale(&orig,&orig,0.25);
                        set_nn_dist(moldyn,0.25*sqrt(3.0)*lc);
-                       ret=diamond_init(a,b,c,lc,atom,&orig);
+                       ret=diamond_init(a,b,c,lc,atom,&orig,p_params,d_params);
                        strcpy(name,"diamond");
                        break;
                default:
@@ -561,25 +613,62 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
 
        /* debug */
        if(ret!=new) {
-               printf("[moldyn] creating lattice failed\n");
+               printf("[moldyn] creating %s lattice (lc=%f) incomplete\n",
+                      name,lc);
+               printf("  (ignore for partial lattice creation)\n");
                printf("  amount of atoms\n");
                printf("  - expected: %d\n",new);
                printf("  - created: %d\n",ret);
+       }
+
+       moldyn->count+=ret;
+       if(ret==new)
+               printf("[moldyn] created %s lattice with %d atoms\n",name,ret);
+
+       for(new=0;new<ret;new++) {
+               atom[new].element=element;
+               atom[new].mass=pse_mass[element];
+               atom[new].attr=attr;
+               atom[new].brand=brand;
+               atom[new].tag=count+new;
+               check_per_bound(moldyn,&(atom[new].r));
+               atom[new].r_0=atom[new].r;
+#ifdef PTHREADS
+               pthread_mutex_init(&(mutex[new]),NULL);
+#endif
+               if(d_params->type) {
+                       new+=1;
+                       atom[new].element=d_params->element;
+                       atom[new].mass=pse_mass[d_params->element];
+                       atom[new].attr=d_params->attr;
+                       atom[new].brand=d_params->brand;
+                       atom[new].tag=count+new;
+                       check_per_bound(moldyn,&(atom[new].r));
+                       atom[new].r_0=atom[new].r;
+#ifdef PTHREADS
+                       pthread_mutex_init(&(mutex[new]),NULL);
+#endif
+               }
+       }
+
+       /* fix allocation */
+       ptr=realloc(moldyn->atom,moldyn->count*sizeof(t_atom));
+       if(!ptr) {
+               perror("[moldyn] realloc (create lattice - alloc fix)");
                return -1;
        }
+       moldyn->atom=ptr;
 
-       moldyn->count+=new;
-       printf("[moldyn] created %s lattice with %d atoms\n",name,new);
+// WHAT ABOUT AMUTEX !!!!
 
-       for(ret=0;ret<new;ret++) {
-               atom[ret].element=element;
-               atom[ret].mass=mass;
-               atom[ret].attr=attr;
-               atom[ret].brand=brand;
-               atom[ret].tag=count+ret;
-               check_per_bound(moldyn,&(atom[ret].r));
-               atom[ret].r_0=atom[ret].r;
+#ifdef LOWMEM_LISTS
+       ptr=realloc(moldyn->lc.subcell->list,moldyn->count*sizeof(int));
+       if(!ptr) {
+               perror("[moldyn] list realloc (create lattice)");
+               return -1;
        }
+       moldyn->lc.subcell->list=ptr;
+#endif
 
        /* update total system mass */
        total_mass_calc(moldyn);
@@ -587,7 +676,7 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        return ret;
 }
 
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
+int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v) {
 
        t_atom *atom;
@@ -613,6 +702,16 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        moldyn->lc.subcell->list=ptr;
 #endif
 
+#ifdef PTHREADS
+       ptr=realloc(amutex,(count+1)*sizeof(pthread_mutex_t));
+       if(!ptr) {
+               perror("[moldyn] mutex realloc (add atom)");
+               return -1;
+       }
+       amutex=ptr;
+       pthread_mutex_init(&(amutex[count]),NULL);
+#endif
+
        atom=moldyn->atom;
 
        /* initialize new atom */
@@ -620,7 +719,7 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
        atom[count].r=*r;
        atom[count].v=*v;
        atom[count].element=element;
-       atom[count].mass=mass;
+       atom[count].mass=pse_mass[element];
        atom[count].brand=brand;
        atom[count].tag=count;
        atom[count].attr=attr;
@@ -662,13 +761,56 @@ int del_atom(t_moldyn *moldyn,int tag) {
        return 0;
 }
 
+#define        set_atom_positions(pos) \
+       if(d_params->type) {\
+               d_o.x=0; d_o.y=0; d_o.z=0;\
+               d_d.x=0; d_d.y=0; d_d.z=0;\
+               switch(d_params->stype) {\
+                       case DEFECT_STYPE_DB_X:\
+                               d_o.x=d_params->od;\
+                               d_d.x=d_params->dd;\
+                               break;\
+                       case DEFECT_STYPE_DB_Y:\
+                               d_o.y=d_params->od;\
+                               d_d.y=d_params->dd;\
+                               break;\
+                       case DEFECT_STYPE_DB_Z:\
+                               d_o.z=d_params->od;\
+                               d_d.z=d_params->dd;\
+                               break;\
+                       case DEFECT_STYPE_DB_R:\
+                               break;\
+                       default:\
+                               printf("[moldyn] WARNING: unknown defect\n");\
+                               break;\
+               }\
+               v3_add(&dr,&pos,&d_o);\
+               v3_copy(&(atom[count].r),&dr);\
+               count+=1;\
+               v3_add(&dr,&pos,&d_d);\
+               v3_copy(&(atom[count].r),&dr);\
+               count+=1;\
+       }\
+       else {\
+               v3_copy(&(atom[count].r),&pos);\
+               count+=1;\
+       }
+
 /* cubic init */
-int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+               t_part_params *p_params,t_defect_params *d_params) {
 
        int count;
        t_3dvec r;
        int i,j,k;
        t_3dvec o;
+       t_3dvec dist;
+       t_3dvec p;
+       t_3dvec d_o;
+       t_3dvec d_d;
+       t_3dvec dr;
+
+       p.x=0; p.y=0; p.z=0;
 
        count=0;
        if(origin)
@@ -676,14 +818,54 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
        else
                v3_zero(&o);
 
+       /* shift partition values */
+       if(p_params->type) {
+               p.x=p_params->p.x+(a*lc)/2.0;
+               p.y=p_params->p.y+(b*lc)/2.0;
+               p.z=p_params->p.z+(c*lc)/2.0;
+       }
+
        r.x=o.x;
        for(i=0;i<a;i++) {
                r.y=o.y;
                for(j=0;j<b;j++) {
                        r.z=o.z;
                        for(k=0;k<c;k++) {
-                               v3_copy(&(atom[count].r),&r);
-                               count+=1;
+                               switch(p_params->type) {
+                                       case PART_INSIDE_R:
+                                               v3_sub(&dist,&r,&p);
+                       if(v3_absolute_square(&dist)<
+                           (p_params->r*p_params->r)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_OUTSIDE_R:
+                                               v3_sub(&dist,&r,&p);
+                       if(v3_absolute_square(&dist)>=
+                           (p_params->r*p_params->r)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)<p_params->d.x)&&
+                          (fabs(dist.y)<p_params->d.y)&&
+                          (fabs(dist.z)<p_params->d.z)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)>=p_params->d.x)||
+                          (fabs(dist.y)>=p_params->d.y)||
+                          (fabs(dist.z)>=p_params->d.z)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       default:        
+                                               set_atom_positions(r);
+                                               break;
+                               }
                                r.z+=lc;
                        }
                        r.y+=lc;
@@ -701,12 +883,18 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 }
 
 /* fcc lattice init */
-int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+             t_part_params *p_params,t_defect_params *d_params) {
 
        int count;
        int i,j,k,l;
        t_3dvec o,r,n;
        t_3dvec basis[3];
+       t_3dvec dist;
+       t_3dvec p;
+       t_3dvec d_d,d_o,dr;
+
+       p.x=0; p.y=0; p.z=0;
 
        count=0;
        if(origin)
@@ -723,6 +911,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
        basis[2].y=0.5*lc;
        basis[2].z=0.5*lc;
 
+       /* shift partition values */
+       if(p_params->type) {
+               p.x=p_params->p.x+(a*lc)/2.0;
+               p.y=p_params->p.y+(b*lc)/2.0;
+               p.z=p_params->p.z+(c*lc)/2.0;
+       }
+
        /* fill up the room */
        r.x=o.x;
        for(i=0;i<a;i++) {
@@ -731,15 +926,81 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
                        r.z=o.z;
                        for(k=0;k<c;k++) {
                                /* first atom */
-                               v3_copy(&(atom[count].r),&r);
-                               count+=1;
-                               r.z+=lc;
+                               switch(p_params->type) {
+                                       case PART_INSIDE_R:
+                                               v3_sub(&dist,&r,&p);
+                       if(v3_absolute_square(&dist)<
+                          (p_params->r*p_params->r)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_OUTSIDE_R:
+                                               v3_sub(&dist,&r,&p);
+                       if(v3_absolute_square(&dist)>=
+                          (p_params->r*p_params->r)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)<p_params->d.x)&&
+                          (fabs(dist.y)<p_params->d.y)&&
+                          (fabs(dist.z)<p_params->d.z)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)>=p_params->d.x)||
+                          (fabs(dist.y)>=p_params->d.y)||
+                          (fabs(dist.z)>=p_params->d.z)) {
+                               set_atom_positions(r);
+                       }
+                                               break;
+                                       default:
+                                               set_atom_positions(r);
+                                               break;
+                               }
                                /* the three face centered atoms */
                                for(l=0;l<3;l++) {
                                        v3_add(&n,&r,&basis[l]);
-                                       v3_copy(&(atom[count].r),&n);
-                                       count+=1;
+                                       switch(p_params->type) {
+                                               case PART_INSIDE_R:
+                       v3_sub(&dist,&n,&p);
+                       if(v3_absolute_square(&dist)<
+                          (p_params->r*p_params->r)) {
+                               set_atom_positions(n);
+                       }
+                                                       break;
+                                               case PART_OUTSIDE_R:
+                       v3_sub(&dist,&n,&p);
+                       if(v3_absolute_square(&dist)>=
+                          (p_params->r*p_params->r)) {
+                               set_atom_positions(n);
+                       }
+                                                       break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&n,&p);
+                       if((fabs(dist.x)<p_params->d.x)&&
+                          (fabs(dist.y)<p_params->d.y)&&
+                          (fabs(dist.z)<p_params->d.z)) {
+                               set_atom_positions(n);
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&n,&p);
+                       if((fabs(dist.x)>=p_params->d.x)||
+                          (fabs(dist.y)>=p_params->d.y)||
+                          (fabs(dist.z)>=p_params->d.z)) {
+                               set_atom_positions(n);
+                       }
+                                               break;
+                                               default:
+                                                       set_atom_positions(n);
+                                                       break;
+                                       }
                                }
+                               r.z+=lc;
                        }
                        r.y+=lc;
                }
@@ -756,12 +1017,13 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
        return count;
 }
 
-int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
+int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
+                 t_part_params *p_params,t_defect_params *d_params) {
 
        int count;
        t_3dvec o;
 
-       count=fcc_init(a,b,c,lc,atom,origin);
+       count=fcc_init(a,b,c,lc,atom,origin,p_params,d_params);
 
        o.x=0.25*lc;
        o.y=0.25*lc;
@@ -769,7 +1031,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
 
        if(origin) v3_add(&o,&o,origin);
 
-       count+=fcc_init(a,b,c,lc,&atom[count],&o);
+       count+=fcc_init(a,b,c,lc,&atom[count],&o,p_params,d_params);
 
        return count;
 }
@@ -847,7 +1109,9 @@ double temperature_calc(t_moldyn *moldyn) {
 
        /* assume up to date kinetic energy, which is 3/2 N k_B T */
 
-       moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+       if(moldyn->count)
+               moldyn->t=(2.0*moldyn->ekin)/(3.0*K_BOLTZMANN*moldyn->count);
+       else moldyn->t=0.0;
 
        return moldyn->t;
 }
@@ -947,8 +1211,8 @@ double virial_sum(t_moldyn *moldyn) {
        }
 
        /* global virial (absolute coordinates) */
-       virial=&(moldyn->gvir);
-       moldyn->gv=virial->xx+virial->yy+virial->zz;
+       //virial=&(moldyn->gvir);
+       //moldyn->gv=virial->xx+virial->yy+virial->zz;
 
        return moldyn->virial;
 }
@@ -969,9 +1233,16 @@ double pressure_calc(t_moldyn *moldyn) {
        moldyn->p=2.0*moldyn->ekin+moldyn->virial;
        moldyn->p/=(3.0*moldyn->volume);
 
+       //moldyn->px=2.0*moldyn->ekinx+moldyn->vir.xx;
+       //moldyn->px/=moldyn->volume;
+       //moldyn->py=2.0*moldyn->ekiny+moldyn->vir.yy;
+       //moldyn->py/=moldyn->volume;
+       //moldyn->pz=2.0*moldyn->ekinz+moldyn->vir.zz;
+       //moldyn->pz/=moldyn->volume;
+
        /* pressure (absolute coordinates) */
-       moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
-       moldyn->gp/=(3.0*moldyn->volume);
+       //moldyn->gp=2.0*moldyn->ekin+moldyn->gv;
+       //moldyn->gp/=(3.0*moldyn->volume);
 
        return moldyn->p;
 }
@@ -996,11 +1267,11 @@ int average_reset(t_moldyn *moldyn) {
 
        /* virial */
        moldyn->virial_sum=0.0;
-       moldyn->gv_sum=0.0;
+       //moldyn->gv_sum=0.0;
 
        /* pressure */
        moldyn->p_sum=0.0;
-       moldyn->gp_sum=0.0;
+       //moldyn->gp_sum=0.0;
        moldyn->tp_sum=0.0;
 
        return 0;
@@ -1038,14 +1309,14 @@ int average_and_fluctuation_calc(t_moldyn *moldyn) {
        /* virial */
        moldyn->virial_sum+=moldyn->virial;
        moldyn->virial_avg=moldyn->virial_sum/denom;
-       moldyn->gv_sum+=moldyn->gv;
-       moldyn->gv_avg=moldyn->gv_sum/denom;
+       //moldyn->gv_sum+=moldyn->gv;
+       //moldyn->gv_avg=moldyn->gv_sum/denom;
 
        /* pressure */
        moldyn->p_sum+=moldyn->p;
        moldyn->p_avg=moldyn->p_sum/denom;
-       moldyn->gp_sum+=moldyn->gp;
-       moldyn->gp_avg=moldyn->gp_sum/denom;
+       //moldyn->gp_sum+=moldyn->gp;
+       //moldyn->gp_avg=moldyn->gp_sum/denom;
        moldyn->tp_sum+=moldyn->tp;
        moldyn->tp_avg=moldyn->tp_sum/denom;
 
@@ -1199,11 +1470,40 @@ int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z) {
        return 0;
 }
 
+int scale_atoms_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+       int i;
+       t_3dvec *r;
+
+       for(i=0;i<moldyn->count;i++) {
+               r=&(moldyn->atom[i].r);
+               r->x*=x;
+               r->y*=y;
+               r->z*=z;
+       }
+
+       return 0;
+}
+
+int scale_dim_ind(t_moldyn *moldyn,double x,double y,double z) {
+
+       t_3dvec *dim;
+
+       dim=&(moldyn->dim);
+
+       dim->x*=x;
+       dim->y*=y;
+       dim->z*=z;
+
+       return 0;
+}
+
 int scale_volume(t_moldyn *moldyn) {
 
        t_3dvec *dim,*vdim;
        double scale;
        t_linkcell *lc;
+       //double sx,sy,sz;
 
        vdim=&(moldyn->vis.dim);
        dim=&(moldyn->dim);
@@ -1218,9 +1518,21 @@ int scale_volume(t_moldyn *moldyn) {
                scale=pow(moldyn->p/moldyn->p_ref,ONE_THIRD);
        }
 
+
+       /*      
+       sx=1.0-(moldyn->p_ref-moldyn->px)*moldyn->p_tc*moldyn->tau;
+       sy=1.0-(moldyn->p_ref-moldyn->py)*moldyn->p_tc*moldyn->tau;
+       sz=1.0-(moldyn->p_ref-moldyn->pz)*moldyn->p_tc*moldyn->tau;
+       sx=pow(sx,ONE_THIRD);
+       sy=pow(sy,ONE_THIRD);
+       sz=pow(sz,ONE_THIRD);
+       */
+
        /* scale the atoms and dimensions */
        scale_atoms(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
        scale_dim(moldyn,SCALE_DIRECT,scale,TRUE,TRUE,TRUE);
+       //scale_atoms_ind(moldyn,sx,sy,sz);
+       //scale_dim_ind(moldyn,sx,sy,sz);
 
        /* visualize dimensions */
        if(vdim->x!=0) {
@@ -1242,6 +1554,9 @@ int scale_volume(t_moldyn *moldyn) {
                lc->x*=scale;
                lc->y*=scale;
                lc->z*=scale;
+               //lc->x*=sx;
+               //lc->y*=sx;
+               //lc->z*=sy;
        }
 
        return 0;
@@ -1255,10 +1570,16 @@ double e_kin_calc(t_moldyn *moldyn) {
 
        atom=moldyn->atom;
        moldyn->ekin=0.0;
+       //moldyn->ekinx=0.0;
+       //moldyn->ekiny=0.0;
+       //moldyn->ekinz=0.0;
 
        for(i=0;i<moldyn->count;i++) {
                atom[i].ekin=0.5*atom[i].mass*v3_absolute_square(&(atom[i].v));
                moldyn->ekin+=atom[i].ekin;
+               //moldyn->ekinx+=0.5*atom[i].mass*atom[i].v.x*atom[i].v.x;
+               //moldyn->ekiny+=0.5*atom[i].mass*atom[i].v.y*atom[i].v.y;
+               //moldyn->ekinz+=0.5*atom[i].mass*atom[i].v.z*atom[i].v.z;
        }
 
        return moldyn->ekin;
@@ -1624,7 +1945,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        struct timeval t1,t2;
        //double tp;
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
        u8 first,change;
        pthread_t io_thread;
        int ret;
@@ -1684,7 +2005,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
        moldyn->debug=0;
 
        /* zero & init moldyn copy */
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
        memset(&md_copy,0,sizeof(t_moldyn));
        atom_copy=malloc(moldyn->count*sizeof(t_atom));
        if(atom_copy==NULL) {
@@ -1693,6 +2014,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
        }
 #endif
 
+#ifdef PTHREADS
+       printf("##################\n");
+       printf("# USING PTHREADS #\n");
+       printf("##################\n");
+#endif
        /* tell the world */
        printf("[moldyn] integration start, go get a coffee ...\n");
 
@@ -1758,7 +2084,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                                dprintf(moldyn->pfd,
                                        "%f %f %f %f %f %f %f\n",moldyn->time,
                                         moldyn->p/BAR,moldyn->p_avg/BAR,
-                                        moldyn->gp/BAR,moldyn->gp_avg/BAR,
+                                        moldyn->p/BAR,moldyn->p_avg/BAR,
                                         moldyn->tp/BAR,moldyn->tp_avg/BAR);
                        }
                }
@@ -1772,7 +2098,11 @@ int moldyn_integrate(t_moldyn *moldyn) {
                if(v) {
                        if(!(moldyn->total_steps%v)) {
                                dprintf(moldyn->vfd,
-                                       "%f %f\n",moldyn->time,moldyn->volume);
+                                       "%f %f %f %f %f\n",moldyn->time,
+                                                          moldyn->volume,
+                                                          moldyn->dim.x,
+                                                          moldyn->dim.y,
+                                                          moldyn->dim.z);
                        }
                }
                if(s) {
@@ -1792,7 +2122,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                }
                if(a) {
                        if(!(moldyn->total_steps%a)) {
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
        /* check whether thread has not terminated yet */
        if(!first) {
                ret=pthread_join(io_thread,NULL);
@@ -1830,7 +2160,7 @@ int moldyn_integrate(t_moldyn *moldyn) {
                        /* get current time */
                        gettimeofday(&t2,NULL);
 
-printf("\rsched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)",
+printf("sched:%d, steps:%d/%d, T:%4.1f/%4.1f P:%4.1f/%4.1f V:%6.1f (%d)\n",
        sched->count,i,moldyn->total_steps,
        moldyn->t,moldyn->t_avg,
        moldyn->p/BAR,moldyn->p_avg/BAR,
@@ -2470,6 +2800,16 @@ int moldyn_read_save_file(t_moldyn *moldyn,char *file) {
                size-=cnt;
        }
 
+#ifdef PTHREADS
+       amutex=malloc(moldyn->count*sizeof(pthread_mutex_t));
+       if(amutex==NULL) {
+               perror("[moldyn] load save file (mutexes)");
+               return -1;
+       }
+       for(cnt=0;cnt<moldyn->count;cnt++)
+               pthread_mutex_init(&(amutex[cnt]),NULL);
+#endif
+
        // hooks etc ...
 
        return 0;
@@ -2974,7 +3314,7 @@ int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
        return 0;
 }
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
 void *visual_atoms(void *ptr) {
 #else
 int visual_atoms(t_moldyn *moldyn) {
@@ -2987,7 +3327,8 @@ int visual_atoms(t_moldyn *moldyn) {
        t_visual *v;
        t_atom *atom;
        t_vb vb;
-#ifdef PTHREADS
+       t_3dvec strain;
+#ifdef VISUAL_THREAD
        t_moldyn *moldyn;
 
        moldyn=ptr;
@@ -3005,7 +3346,7 @@ int visual_atoms(t_moldyn *moldyn) {
        vb.fd=open(file,O_WRONLY|O_CREAT|O_TRUNC,S_IRUSR|S_IWUSR);
        if(vb.fd<0) {
                perror("open visual save file fd");
-#ifndef PTHREADS
+#ifndef VISUAL_THREAD
                return -1;
 #endif
        }
@@ -3017,17 +3358,23 @@ int visual_atoms(t_moldyn *moldyn) {
                moldyn->count,moldyn->time,help/40.0,help/40.0,-0.8*help);
 
        // atomic configuration
-       for(i=0;i<moldyn->count;i++)
+       for(i=0;i<moldyn->count;i++) {
+               v3_sub(&strain,&(atom[i].r),&(atom[i].r_0));
+               check_per_bound(moldyn,&strain);
                // atom type, positions, color and kinetic energy
                dprintf(vb.fd,"%s %f %f %f %s %f\n",pse_name[atom[i].element],
                                                    atom[i].r.x,
                                                    atom[i].r.y,
                                                    atom[i].r.z,
                                                    pse_col[atom[i].element],
-                                                   atom[i].ekin);
+                                                   //atom[i].ekin);
+                                                   sqrt(v3_absolute_square(&strain)));
+       }
        
        // bonds between atoms
+#ifndef VISUAL_THREAD
        process_2b_bonds(moldyn,&vb,visual_bonds_process);
+#endif
        
        // boundaries
        if(dim.x) {
@@ -3073,7 +3420,7 @@ int visual_atoms(t_moldyn *moldyn) {
 
        close(vb.fd);
 
-#ifdef PTHREADS
+#ifdef VISUAL_THREAD
        pthread_exit(NULL);
 
 }