security checkin: fixes to lattice stage, px py pz calc (now default!

TAKE CARE) and stuff i forgot ...

diff --git a/mdrun.c b/mdrun.c
index 83f501ce2074cb568b7637bc0f8f2c4d2a00146e..6effe1aff6c56f33b635da11e332e27a4252b65e 100644 (file)
--- a/mdrun.c
+++ b/mdrun.c
@@ -122,6 +122,12 @@ int add_stage(t_mdrun *mdrun,u8 type,void *params) {
                case STAGE_SET_TIMESTEP:
                        psize=sizeof(t_set_timestep_params);
                        break;
+               case STAGE_FILL:
+                       psize=sizeof(t_fill_params);
+                       break;
+               case STAGE_THERMAL_INIT:
+                       psize=0;
+                       break;
                default:
                        printf("%s unknown stage type: %02x\n",ME,type);
                        return -1;
@@ -168,6 +174,7 @@ int mdrun_parse_config(t_mdrun *mdrun) {
        t_chsattr_params csp;
        t_set_temp_params stp;
        t_set_timestep_params stsp;
+       t_fill_params fp;
 
        /* open config file */
        fd=open(mdrun->cfile,O_RDONLY);
@@ -199,6 +206,7 @@ int mdrun_parse_config(t_mdrun *mdrun) {
                memset(&csp,0,sizeof(t_chsattr_params));
                memset(&stp,0,sizeof(t_set_temp_params));
                memset(&stsp,0,sizeof(t_set_timestep_params));
+               memset(&fp,0,sizeof(t_fill_params));
 
                // get command + args
                wcnt=0;
@@ -279,19 +287,60 @@ int mdrun_parse_config(t_mdrun *mdrun) {
                        mdrun->m2=pse_mass[mdrun->element2];
                }
                else if(!strncmp(word[0],"fill",6)) {
-                       // only lc mode by now
-                       mdrun->lx=atoi(word[2]);
-                       mdrun->ly=atoi(word[3]);
-                       mdrun->lz=atoi(word[4]);
-                       mdrun->lc=atof(word[5]);
-                       if(wcnt==8) {
-                               mdrun->fill_element=atoi(word[6]);
-                               mdrun->fill_brand=atoi(word[7]);
+                       fp.lx=atoi(word[2]);
+                       fp.ly=atoi(word[3]);
+                       fp.lz=atoi(word[4]);
+                       fp.lc=atof(word[5]);
+                       mdrun->lc=fp.lc;
+                       if(!strncmp(word[1],"lc",2)) {
+                               if(wcnt==8) {
+                                       fp.fill_element=atoi(word[6]);
+                                       fp.fill_brand=atoi(word[7]);
+                               }
+                               else {
+                                       fp.fill_element=mdrun->element1;
+                                       fp.fill_brand=0;
+                               }
                        }
                        else {
-                               mdrun->fill_element=mdrun->element1;
-                               mdrun->fill_brand=0;
+                               switch(word[6][0]) {
+                                       case 'i':
+                                               if(word[6][1]=='r')
+                                                       fp.p_type=PART_INSIDE_R;
+                                               else
+                                                       fp.p_type=PART_INSIDE_D;
+                                               break;
+                                       case 'o':
+                                               if(word[6][1]=='r')
+                                               fp.p_type=PART_OUTSIDE_R;
+                                               else
+                                               fp.p_type=PART_OUTSIDE_D;
+                                               break;
+                                       default:
+                                               break;
+                               }
+                       }
+                       if((fp.p_type==PART_INSIDE_R)||
+                           (fp.p_type==PART_OUTSIDE_R)) {
+                               fp.p_vals.r=atof(word[7]);      
+                               fp.p_vals.p.x=atof(word[8]);    
+                               fp.p_vals.p.y=atof(word[9]);    
+                               fp.p_vals.p.z=atof(word[10]);
                        }
+                       if((fp.p_type==PART_INSIDE_D)||
+                           (fp.p_type==PART_OUTSIDE_D)) {
+                               fp.p_vals.p.x=atof(word[7]);    
+                               fp.p_vals.p.y=atof(word[8]);    
+                               fp.p_vals.p.z=atof(word[9]);
+                               fp.p_vals.d.x=atof(word[10]);   
+                               fp.p_vals.d.y=atof(word[11]);   
+                               fp.p_vals.d.z=atof(word[12]);
+                       }
+                       fp.lattice=mdrun->lattice;
+                       add_stage(mdrun,STAGE_FILL,&fp);
+               }
+               else if(!strncmp(word[0],"thermal_init",12)) {
+                       add_stage(mdrun,STAGE_THERMAL_INIT,NULL);
                }
                else if(!strncmp(word[0],"aattr",5)) {
                        // for aatrib line we need a special stage
@@ -937,6 +986,7 @@ int mdrun_hook(void *ptr1,void *ptr2) {
        t_list *sl;
        int steps;
        u8 change_stage;
+       t_3dvec o;
 
        t_insert_atoms_params *iap;
        t_insert_mixed_atoms_params *imp;
@@ -944,6 +994,7 @@ int mdrun_hook(void *ptr1,void *ptr2) {
        t_anneal_params *ap;
        t_set_temp_params *stp;
        t_set_timestep_params *stsp;
+       t_fill_params *fp;
 
        moldyn=ptr1;
        mdrun=ptr2;
@@ -1049,6 +1100,51 @@ int mdrun_hook(void *ptr1,void *ptr2) {
                                mdrun->timestep=stsp->tau;
                                change_stage=TRUE;
                                break;
+                       case STAGE_FILL:
+                               stage_print("  -> fill lattice\n\n");
+                               fp=stage->params;
+                               if(fp->lattice!=ZINCBLENDE) {
+                                       create_lattice(moldyn,
+                                                      fp->lattice,fp->lc,
+                                                      mdrun->element1,
+                                                      mdrun->m1,
+                                                      DEFAULT_ATOM_ATTR,0,
+                                                      fp->lx,fp->ly,fp->lz,
+                                                      NULL,fp->p_type,
+                                                      &(fp->p_vals));
+                               }
+                               else {
+                                       o.x=0.5*0.25*fp->lc;
+                                       o.y=o.x;
+                                       o.z=o.x;
+                                       create_lattice(moldyn,
+                                                      FCC,fp->lc,
+                                                      mdrun->element1,
+                                                      mdrun->m1,
+                                                      DEFAULT_ATOM_ATTR,0,
+                                                      fp->lx,fp->ly,fp->lz,
+                                                      &o,fp->p_type,
+                                                      &(fp->p_vals));
+                                       o.x+=0.25*fp->lc;
+                                       o.y=o.x;
+                                       o.z=o.x;
+                                       create_lattice(moldyn,
+                                                      FCC,fp->lc,
+                                                      mdrun->element2,
+                                                      mdrun->m2,
+                                                      DEFAULT_ATOM_ATTR,1,
+                                                      fp->lx,fp->ly,fp->lz,
+                                                      &o,fp->p_type,
+                                                      &(fp->p_vals));
+                               }
+                               moldyn_bc_check(moldyn);
+                               change_stage=TRUE;
+                               break;
+                       case STAGE_THERMAL_INIT:
+                               stage_print("  -> thermal init\n\n");
+                               thermal_init(moldyn,TRUE);
+                               change_stage=TRUE;
+                               break;
                        default:
                                printf("%s unknwon stage type\n",ME);
                                break;
@@ -1086,7 +1182,7 @@ int main(int argc,char **argv) {
 
        t_mdrun mdrun;
        t_moldyn moldyn;
-       t_3dvec o;
+       //t_3dvec o;
 
        /* clear structs */
        memset(&mdrun,0,sizeof(t_mdrun));
@@ -1161,6 +1257,7 @@ int main(int argc,char **argv) {
        /* initial lattice and dimensions */
        set_dim(&moldyn,mdrun.dim.x,mdrun.dim.y,mdrun.dim.z,mdrun.vis);
        set_pbc(&moldyn,mdrun.pbcx,mdrun.pbcy,mdrun.pbcz);
+       /* replaced by fill stage !! 
        switch(mdrun.lattice) {
                case FCC:
                        create_lattice(&moldyn,FCC,mdrun.lc,mdrun.fill_element,
@@ -1194,11 +1291,14 @@ int main(int argc,char **argv) {
                        return -1;
        }
        moldyn_bc_check(&moldyn);
+       replaced by fill stage !! */
 
        /* temperature and pressure */
        set_temperature(&moldyn,mdrun.temperature);
        set_pressure(&moldyn,mdrun.pressure);
+       /* replaced thermal_init stage
        thermal_init(&moldyn,TRUE);
+       */
 
 addsched:
        /* first schedule */

diff --git a/mdrun.h b/mdrun.h
index 71c4d83d36ef95e67b17ef7e5d6c7c89349427b2..d6476e541e5db74d14dbdc70766510fe0eb517fd 100644 (file)
--- a/mdrun.h
+++ b/mdrun.h
@@ -50,6 +50,8 @@ typedef struct s_stage {
 #define STAGE_CHSATTR                          0x06
 #define STAGE_SET_TEMP                         0x07
 #define STAGE_SET_TIMESTEP                     0x08
+#define STAGE_FILL                             0x09
+#define STAGE_THERMAL_INIT                     0x10
 
 typedef struct s_mdrun {
        char cfile[128];                        // config file
@@ -72,15 +74,10 @@ typedef struct s_mdrun {
        double m1;
        int element2;                           // element 2
        double m2;
+
        double lc;                              // lattice constant
-       int lx;                                 // amount of lc units
-       int ly;
-       int lz;
        u8 lattice;                             // type of lattice
 
-       int fill_element;
-       u8 fill_brand;
-
        u8 sattr;                               // system attributes
        double temperature;                     // temperature
        double pressure;                        // pressure
@@ -194,9 +191,20 @@ typedef struct s_set_timestep_params {
        double tau;
 } t_set_timestep_params;
 
+typedef struct s_fill_params {
+       double lc;                              // lattice constant
+       int lx;                                 // amount of lc units
+       int ly;
+       int lz;
+       u8 lattice;
+       int fill_element;
+       u8 fill_brand;
+       u8 p_type;
+       t_part_vals p_vals;
+} t_fill_params;
+
 /*
  * function prototypes
  */
 
-
 #endif

diff --git a/moldyn.c b/moldyn.c
index 8451c4436f43db2213bec83677c359fb4062f133..a9341941909c452a613ff62877e6b016bbf10571 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -596,29 +596,48 @@ 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 in case of partial lattice creation)\n");
                printf("  amount of atoms\n");
                printf("  - expected: %d\n",new);
                printf("  - created: %d\n",ret);
-               return -1;
        }
 
-       moldyn->count+=new;
-       printf("[moldyn] created %s lattice with %d atoms\n",name,new);
+       moldyn->count+=ret;
+       if(ret==new)
+               printf("[moldyn] created %s lattice with %d atoms\n",name,ret);
 
-       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;
+       for(new=0;new<ret;new++) {
+               atom[new].element=element;
+               atom[new].mass=mass;
+               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[ret]),NULL);
+               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;
+
+#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);
 
@@ -719,6 +738,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
        int i,j,k;
        t_3dvec o;
        t_3dvec dist;
+       t_3dvec p;
+
+       p.x=0; p.y=0; p.z=0;
 
        count=0;
        if(origin)
@@ -728,9 +750,9 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
        /* shift partition values */
        if(p_type) {
-               p_vals->p.x+=(a*lc)/2.0;
-               p_vals->p.y+=(b*lc)/2.0;
-               p_vals->p.z+=(c*lc)/2.0;
+               p.x=p_vals->p.x+(a*lc)/2.0;
+               p.y=p_vals->p.y+(b*lc)/2.0;
+               p.z=p_vals->p.z+(c*lc)/2.0;
        }
 
        r.x=o.x;
@@ -741,17 +763,35 @@ int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
                        for(k=0;k<c;k++) {
                                switch(p_type) {
                                        case PART_INSIDE_R:
-                                               v3_sub(&dist,&r,&(p_vals->p));
+                                               v3_sub(&dist,&r,&p);
                        if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
                                v3_copy(&(atom[count].r),&r);
                                count+=1;
                        }
                                                break;
                                        case PART_OUTSIDE_R:
-                                               v3_sub(&dist,&r,&(p_vals->p));
+                                               v3_sub(&dist,&r,&p);
                        if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
                                v3_copy(&(atom[count].r),&r);
                                count+=1;
+                       }
+                                               break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)<p_vals->d.x)&&
+                          (fabs(dist.y)<p_vals->d.y)&&
+                          (fabs(dist.z)<p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)>=p_vals->d.x)||
+                          (fabs(dist.y)>=p_vals->d.y)||
+                          (fabs(dist.z)>=p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
                        }
                                                break;
                                        default:        
@@ -784,6 +824,9 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
        t_3dvec o,r,n;
        t_3dvec basis[3];
        t_3dvec dist;
+       t_3dvec p;
+
+       p.x=0; p.y=0; p.z=0;
 
        count=0;
        if(origin)
@@ -802,9 +845,9 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
 
        /* shift partition values */
        if(p_type) {
-               p_vals->p.x+=(a*lc)/2.0;
-               p_vals->p.y+=(b*lc)/2.0;
-               p_vals->p.z+=(c*lc)/2.0;
+               p.x=p_vals->p.x+(a*lc)/2.0;
+               p.y=p_vals->p.y+(b*lc)/2.0;
+               p.z=p_vals->p.z+(c*lc)/2.0;
        }
 
        /* fill up the room */
@@ -817,17 +860,35 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
                                /* first atom */
                                switch(p_type) {
                                        case PART_INSIDE_R:
-                                               v3_sub(&dist,&r,&(p_vals->p));
+                                               v3_sub(&dist,&r,&p);
                        if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
                                v3_copy(&(atom[count].r),&r);
                                count+=1;
                        }
                                                break;
                                        case PART_OUTSIDE_R:
-                                               v3_sub(&dist,&r,&(p_vals->p));
+                                               v3_sub(&dist,&r,&p);
                        if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
                                v3_copy(&(atom[count].r),&r);
                                count+=1;
+                       }
+                                               break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)<p_vals->d.x)&&
+                          (fabs(dist.y)<p_vals->d.y)&&
+                          (fabs(dist.z)<p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&r,&p);
+                       if((fabs(dist.x)>=p_vals->d.x)||
+                          (fabs(dist.y)>=p_vals->d.y)||
+                          (fabs(dist.z)>=p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
                        }
                                                break;
                                        default:
@@ -840,19 +901,37 @@ int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin,
                                        v3_add(&n,&r,&basis[l]);
                                        switch(p_type) {
                                                case PART_INSIDE_R:
-                       v3_sub(&dist,&n,&(p_vals->p));
+                       v3_sub(&dist,&n,&p);
                        if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
-                               v3_copy(&(atom[count].r),&r);
+                               v3_copy(&(atom[count].r),&n);
                                count+=1;
                        }
                                                        break;
                                                case PART_OUTSIDE_R:
-                       v3_sub(&dist,&n,&(p_vals->p));
+                       v3_sub(&dist,&n,&p);
                        if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
-                               v3_copy(&(atom[count].r),&r);
+                               v3_copy(&(atom[count].r),&n);
                                count+=1;
                        }
                                                        break;
+                                       case PART_INSIDE_D:
+                                               v3_sub(&dist,&n,&p);
+                       if((fabs(dist.x)<p_vals->d.x)&&
+                          (fabs(dist.y)<p_vals->d.y)&&
+                          (fabs(dist.z)<p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&n);
+                               count+=1;
+                       }
+                                               break;
+                                       case PART_OUTSIDE_D:
+                                               v3_sub(&dist,&n,&p);
+                       if((fabs(dist.x)>=p_vals->d.x)||
+                          (fabs(dist.y)>=p_vals->d.y)||
+                          (fabs(dist.z)>=p_vals->d.z)) {
+                               v3_copy(&(atom[count].r),&n);
+                               count+=1;
+                       }
+                                               break;
                                                default:
                                        v3_copy(&(atom[count].r),&n);
                                        count+=1;
@@ -968,7 +1047,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;
 }
@@ -1068,8 +1149,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;
 }
@@ -1090,9 +1171,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;
 }
@@ -1117,11 +1205,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;
@@ -1159,14 +1247,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;
 
@@ -1320,17 +1408,47 @@ 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;
+       //double scale;
        t_linkcell *lc;
+       double sx,sy,sz;
 
        vdim=&(moldyn->vis.dim);
        dim=&(moldyn->dim);
        lc=&(moldyn->lc);
 
        /* scaling factor */
+       /*
        if(moldyn->pt_scale&P_SCALE_BERENDSEN) {
                scale=1.0-(moldyn->p_ref-moldyn->p)*moldyn->p_tc*moldyn->tau;
                scale=pow(scale,ONE_THIRD);
@@ -1338,10 +1456,20 @@ int scale_volume(t_moldyn *moldyn) {
        else {
                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(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) {
@@ -1360,9 +1488,12 @@ int scale_volume(t_moldyn *moldyn) {
                link_cell_shutdown(moldyn);
                link_cell_init(moldyn,QUIET);
        } else {
-               lc->x*=scale;
-               lc->y*=scale;
-               lc->z*=scale;
+               //lc->x*=scale;
+               //lc->y*=scale;
+               //lc->z*=scale;
+               lc->x*=sx;
+               lc->y*=sx;
+               lc->z*=sy;
        }
 
        return 0;
@@ -1376,10 +1507,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;
@@ -1884,7 +2021,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);
                        }
                }
@@ -1898,7 +2035,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) {
@@ -2596,6 +2737,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;

diff --git a/moldyn.h b/moldyn.h
index 1e9db2a3ae3527e17b33b37dd03b641c5ddceeb9..d065d43e2d3420a79df50ae334f68a8ff43800dc 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -137,14 +137,22 @@ typedef struct s_moldyn {
        double t_sum;           /* sum over all t */
        double t_avg;           /* average value of t */
 
+       /* for sale! */
        t_virial gvir;          /* global virial (absolute coordinates) */
-       double gv;
-       double gv_sum;
-       double gv_avg;
-
-       double gp;              /* pressure computed from global virial */
-       double gp_sum;          /* sum over all gp */
-       double gp_avg;          /* average value of gp */
+       //double gv;
+       //double gv_sum;
+       //double gv_avg;
+       double sale1;
+       double sale2;
+       double sale3;
+
+       // gp stuff exchanged by kinetic energies
+       //double gp;            /* pressure computed from global virial */
+       //double gp_sum;                /* sum over all gp */
+       //double gp_avg;                /* average value of gp */
+       double ekinx;
+       double ekiny;
+       double ekinz;
 
        t_virial vir;           /* actual virial */
        double virial;
@@ -258,6 +266,8 @@ typedef struct s_part_vals {
 
 #define PART_INSIDE_R   1
 #define PART_OUTSIDE_R  2
+#define PART_INSIDE_D   3
+#define PART_OUTSIDE_D  4
 
 /*
  *