actually just bullshit ... (added unique tag -> maybe symmetry)

diff --git a/moldyn.c b/moldyn.c
index 551e07a..f8bfd3f 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -230,33 +230,38 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  */
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 bnum,int a,int b,int c) {
+                   u8 attr,u8 brand,int a,int b,int c) {
 
-       int count;
+       int new,count;
        int ret;
        t_3dvec origin;
+       void *ptr;
+       t_atom *atom;
 
-       count=a*b*c;
+       new=a*b*c;
+       count=moldyn->count;
 
        /* how many atoms do we expect */
-       if(type==FCC) count*=4;
-       if(type==DIAMOND) count*=8;
+       if(type==FCC) new*=4;
+       if(type==DIAMOND) new*=8;
 
        /* allocate space for atoms */
-       moldyn->atom=malloc(count*sizeof(t_atom));
-       if(moldyn->atom==NULL) {
-               perror("malloc (atoms)");
+       ptr=realloc(moldyn->atom,(count+new)*sizeof(t_atom));
+       if(!ptr) {
+               perror("[moldyn] realloc (create lattice)");
                return -1;
        }
-
+       moldyn->atom=ptr;
+       atom=&(moldyn->atom[count]);
+               
        v3_zero(&origin);
 
        switch(type) {
                case FCC:
-                       ret=fcc_init(a,b,c,lc,moldyn->atom,&origin);
+                       ret=fcc_init(a,b,c,lc,atom,&origin);
                        break;
                case DIAMOND:
-                       ret=diamond_init(a,b,c,lc,moldyn->atom,&origin);
+                       ret=diamond_init(a,b,c,lc,atom,&origin);
                        break;
                default:
                        printf("unknown lattice type (%02x)\n",type);
@@ -264,24 +269,24 @@ int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
        }
 
        /* debug */
-       if(ret!=count) {
+       if(ret!=new) {
                printf("[moldyn] creating lattice failed\n");
                printf("  amount of atoms\n");
-               printf("  - expected: %d\n",count);
+               printf("  - expected: %d\n",new);
                printf("  - created: %d\n",ret);
                return -1;
        }
 
-       moldyn->count=count;
-       printf("[moldyn] created lattice with %d atoms\n",count);
+       moldyn->count+=new;
+       printf("[moldyn] created lattice with %d atoms\n",new);
 
-       while(count) {
-               count-=1;
-               moldyn->atom[count].element=element;
-               moldyn->atom[count].mass=mass;
-               moldyn->atom[count].attr=attr;
-               moldyn->atom[count].bnum=bnum;
-               check_per_bound(moldyn,&(moldyn->atom[count].r));
+       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));
        }
 
        return ret;
@@ -371,7 +376,7 @@ int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
        return count;
 }
 
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v) {
 
        t_atom *atom;
@@ -379,9 +384,9 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
        int count;
        
        atom=moldyn->atom;
-       count=++(moldyn->count);
+       count=(moldyn->count)++;
 
-       ptr=realloc(atom,count*sizeof(t_atom));
+       ptr=realloc(atom,(count+1)*sizeof(t_atom));
        if(!ptr) {
                perror("[moldyn] realloc (add atom)");
                return -1;
@@ -389,12 +394,13 @@ int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
        moldyn->atom=ptr;
 
        atom=moldyn->atom;
-       atom[count-1].r=*r;
-       atom[count-1].v=*v;
-       atom[count-1].element=element;
-       atom[count-1].mass=mass;
-       atom[count-1].bnum=bnum;
-       atom[count-1].attr=attr;
+       atom[count].r=*r;
+       atom[count].v=*v;
+       atom[count].element=element;
+       atom[count].mass=mass;
+       atom[count].brand=brand;
+       atom[count].tag=count;
+       atom[count].attr=attr;
 
        return 0;
 }
@@ -1243,11 +1249,11 @@ int tersoff_mult_complete_params(t_tersoff_mult_params *p) {
 /* tersoff 1 body part */
 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
 
-       int num;
+       int brand;
        t_tersoff_mult_params *params;
        t_tersoff_exchange *exchange;
        
-       num=ai->bnum;
+       brand=ai->brand;
        params=moldyn->pot1b_params;
        exchange=&(params->exchange);
 
@@ -1256,15 +1262,15 @@ int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai) {
         *         their right values
         */
 
-       exchange->beta_i=&(params->beta[num]);
-       exchange->n_i=&(params->n[num]);
-       exchange->c_i=&(params->c[num]);
-       exchange->d_i=&(params->d[num]);
-       exchange->h_i=&(params->h[num]);
+       exchange->beta_i=&(params->beta[brand]);
+       exchange->n_i=&(params->n[brand]);
+       exchange->c_i=&(params->c[brand]);
+       exchange->d_i=&(params->d[brand]);
+       exchange->h_i=&(params->h[brand]);
 
        exchange->betaini=pow(*(exchange->beta_i),*(exchange->n_i));
-       exchange->ci2=params->c[num]*params->c[num];
-       exchange->di2=params->d[num]*params->d[num];
+       exchange->ci2=params->c[brand]*params->c[brand];
+       exchange->di2=params->d[brand]*params->d[brand];
        exchange->ci2di2=exchange->ci2/exchange->di2;
 
        return 0;
@@ -1280,12 +1286,12 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        double A,B,R,S,lambda,mu;
        double f_r,df_r;
        double f_c,df_c;
-       int num;
+       int brand;
        double s_r;
        double arg;
 
        params=moldyn->pot2b_params;
-       num=aj->bnum;
+       brand=aj->brand;
        exchange=&(params->exchange);
 
        /* clear 3bp and 2bp post run */
@@ -1318,13 +1324,13 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        exchange->dist_ij=dist_ij;
 
        /* constants */
-       if(num==ai->bnum) {
-               S=params->S[num];
-               R=params->R[num];
-               A=params->A[num];
-               B=params->B[num];
-               lambda=params->lambda[num];
-               mu=params->mu[num];
+       if(brand==ai->brand) {
+               S=params->S[brand];
+               R=params->R[brand];
+               A=params->A[brand];
+               B=params->B[brand];
+               lambda=params->lambda[brand];
+               mu=params->mu[brand];
                exchange->chi=1.0;
        }
        else {
@@ -1342,12 +1348,12 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
                return 0;
 
        /* more constants */
-       exchange->beta_j=&(params->beta[num]);
-       exchange->n_j=&(params->n[num]);
-       exchange->c_j=&(params->c[num]);
-       exchange->d_j=&(params->d[num]);
-       exchange->h_j=&(params->h[num]);
-       if(num==ai->bnum) {
+       exchange->beta_j=&(params->beta[brand]);
+       exchange->n_j=&(params->n[brand]);
+       exchange->c_j=&(params->c[brand]);
+       exchange->d_j=&(params->d[brand]);
+       exchange->h_j=&(params->h[brand]);
+       if(brand==ai->brand) {
                exchange->betajnj=exchange->betaini;
                exchange->cj2=exchange->ci2;
                exchange->dj2=exchange->di2;
@@ -1355,8 +1361,8 @@ int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) {
        }
        else {
                exchange->betajnj=pow(*(exchange->beta_j),*(exchange->n_j));
-               exchange->cj2=params->c[num]*params->c[num];
-               exchange->dj2=params->d[num]*params->d[num];
+               exchange->cj2=params->c[brand]*params->c[brand];
+               exchange->dj2=params->d[brand]*params->d[brand];
                exchange->cj2dj2=exchange->cj2/exchange->dj2;
        }
 
@@ -1536,7 +1542,7 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        double h_cos,d2_h_cos2;
        double frac,g,zeta,chi;
        double tmp;
-       int num;
+       int brand;
 
        params=moldyn->pot3b_params;
        exchange=&(params->exchange);
@@ -1588,10 +1594,10 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d_ik=v3_norm(&dist_ik);
 
        /* ik constants */
-       num=ai->bnum;
-       if(num==ak->bnum) {
-               R=params->R[num];
-               S=params->S[num];
+       brand=ai->brand;
+       if(brand==ak->brand) {
+               R=params->R[brand];
+               S=params->S[brand];
        }
        else {
                R=params->Rmixed;
@@ -1674,12 +1680,12 @@ int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc) {
        d_jk=v3_norm(&dist_jk);
 
        /* jk constants */
-       num=aj->bnum;
-       if(num==ak->bnum) {
-               R=params->R[num];
-               S=params->S[num];
-               B=params->B[num];
-               mu=params->mu[num];
+       brand=aj->brand;
+       if(brand==ak->brand) {
+               R=params->R[brand];
+               S=params->S[brand];
+               B=params->B[brand];
+               mu=params->mu[brand];
                chi=1.0;
        }
        else {
@@ -1787,13 +1793,13 @@ int moldyn_bc_check(t_moldyn *moldyn) {
        t_atom *atom;
        t_3dvec *dim;
        int i;
-double x;
-u8 byte;
-int j,k;
+       double x;
+       u8 byte;
+       int j,k;
 
        atom=moldyn->atom;
        dim=&(moldyn->dim);
-x=dim->x/2;
+       x=dim->x/2;
 
        for(i=0;i<moldyn->count;i++) {
                if(atom[i].r.x>=dim->x/2||-atom[i].r.x>dim->x/2) {

diff --git a/moldyn.h b/moldyn.h
index b4f71d6..9493f5f 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -27,11 +27,12 @@ typedef struct s_atom {
        t_3dvec r;              /* position */
        t_3dvec v;              /* velocity */
        t_3dvec f;              /* force */
-       t_3dvec virial;         /* virial */
+       t_3dvec virial;         /* virial (v_xx, v_yy, v_zz) */
        double e;               /* site energy */
        int element;            /* number of element in pse */
        double mass;            /* atom mass */
-       u8 bnum;                /* brand number */
+       u8 brand;               /* brand id */
+       int tag;                /* atom unique id (number of atom) */
        u8 attr;                /* attributes */
 } t_atom;
 
@@ -377,10 +378,10 @@ int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
 int moldyn_log_shutdown(t_moldyn *moldyn);
 
 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
-                   u8 attr,u8 bnum,int a,int b,int c);
+                   u8 attr,u8 brand,int a,int b,int c);
 int fcc_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);
-int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
+int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr,
              t_3dvec *r,t_3dvec *v);
 int destroy_atoms(t_moldyn *moldyn);
 

diff --git a/sic.c b/sic.c
index a105987..f4cdd7d 100644 (file)
--- a/sic.c
+++ b/sic.c
@@ -137,8 +137,8 @@ int main(int argc,char **argv) {
 
        /* set temperature */
        printf("[sic] setting temperature\n");
-       set_temperature(&md,273.0+1410.0);
-       //set_temperature(&md,273.0+450.0);
+       //set_temperature(&md,273.0+1410.0);
+       set_temperature(&md,273.0+450.0);
        //set_temperature(&md,273.0);
        //set_temperature(&md,1.0);
        //set_temperature(&md,0.0);
@@ -151,6 +151,7 @@ int main(int argc,char **argv) {
        printf("[sic] set p/t scaling\n");
        //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
        //                 T_SCALE_BERENDSEN,100.0);
+       set_pt_scale(&md,0,0,T_SCALE_BERENDSEN,100.0);
        
        /* initial thermal fluctuations of particles (in equilibrium) */
        printf("[sic] thermal init\n");