added partition support in create_lattice function

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 13d7e5d..8451c44 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -513,7 +513,8 @@ int moldyn_log_shutdown(t_moldyn *moldyn) {
  */
 
 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) {
+                   u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin,
+                   u8 p_type,t_part_vals *p_vals) {
 
        int new,count;
        int ret;
@@ -571,21 +572,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_type,p_vals);
                        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_type,p_vals);
                        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_type,p_vals);
                        strcpy(name,"diamond");
                        break;
                default:
@@ -710,12 +711,14 @@ int del_atom(t_moldyn *moldyn,int tag) {
 }
 
 /* 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,
+               u8 p_type,t_part_vals *p_vals) {
 
        int count;
        t_3dvec r;
        int i,j,k;
        t_3dvec o;
+       t_3dvec dist;
 
        count=0;
        if(origin)
@@ -723,14 +726,39 @@ 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_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;
+       }
+
        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++) {
+                               switch(p_type) {
+                                       case PART_INSIDE_R:
+                                               v3_sub(&dist,&r,&(p_vals->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));
+                       if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
+                       }
+                                               break;
+                                       default:        
+                                               v3_copy(&(atom[count].r),&r);
+                                               count+=1;
+                                               break;
+                               }
                                r.z+=lc;
                        }
                        r.y+=lc;
@@ -748,12 +776,14 @@ 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,
+             u8 p_type,t_part_vals *p_vals) {
 
        int count;
        int i,j,k,l;
        t_3dvec o,r,n;
        t_3dvec basis[3];
+       t_3dvec dist;
 
        count=0;
        if(origin)
@@ -770,6 +800,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_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;
+       }
+
        /* fill up the room */
        r.x=o.x;
        for(i=0;i<a;i++) {
@@ -778,15 +815,51 @@ 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 */
+                               switch(p_type) {
+                                       case PART_INSIDE_R:
+                                               v3_sub(&dist,&r,&(p_vals->p));
+                       if(v3_absolute_square(&dist)<(p_vals->r*p_vals->r)) {
                                v3_copy(&(atom[count].r),&r);
                                count+=1;
-                               r.z+=lc;
+                       }
+                                               break;
+                                       case PART_OUTSIDE_R:
+                                               v3_sub(&dist,&r,&(p_vals->p));
+                       if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
+                       }
+                                               break;
+                                       default:
+                                               v3_copy(&(atom[count].r),&r);
+                                               count+=1;
+                                               break;
+                               }
                                /* the three face centered atoms */
                                for(l=0;l<3;l++) {
                                        v3_add(&n,&r,&basis[l]);
+                                       switch(p_type) {
+                                               case PART_INSIDE_R:
+                       v3_sub(&dist,&n,&(p_vals->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,&n,&(p_vals->p));
+                       if(v3_absolute_square(&dist)>=(p_vals->r*p_vals->r)) {
+                               v3_copy(&(atom[count].r),&r);
+                               count+=1;
+                       }
+                                                       break;
+                                               default:
                                        v3_copy(&(atom[count].r),&n);
                                        count+=1;
+                                                       break;
+                                       }
                                }
+                               r.z+=lc;
                        }
                        r.y+=lc;
                }
@@ -803,12 +876,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,
+                 u8 p_type,t_part_vals *p_vals) {
 
        int count;
        t_3dvec o;
 
-       count=fcc_init(a,b,c,lc,atom,origin);
+       count=fcc_init(a,b,c,lc,atom,origin,p_type,p_vals);
 
        o.x=0.25*lc;
        o.y=0.25*lc;
@@ -816,7 +890,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_type,p_vals);
 
        return count;
 }