well, sure, just increment each basis by 1 ... (should work now!)
[physik/posic.git] / vasp_tools / create_lattice.c
index 5f6cb8d..4b258a0 100644 (file)
 #include <math.h>
 #include "../math/math.h"
 
-#define FCC            0x01
-#define DIAMOND                0x02
-#define ZINCBLENDE     0x03
+/*
+ * lattice types:
+ *
+ * there is fcc/diamond/zincblende only!
+ *
+ * basis:
+ *
+ * 0: <0.5 0.5 0> <0 0.5 0.5> <0.5 0 0.5> contains 2 atoms
+ * 1: <0.5 -0.5 0> <0.5 0.5 0> <0 0 1> contains 4 atoms
+ * 2: <1 0 0> <0 1 0> <0 0 1> contains 8 atoms
+ *
+ * atoms:
+ *
+ * - type 0:
+ *   1st fcc: <0 0 0>
+ *   2nd fcc: <0.25 0.25 0.25>
+ *   -> diag: <0.25 0.25 0.25>
+ *
+ * - type 1:
+ *   1st fcc: <0 0 0> <0.5 0.5 0.5>
+ *   2nd fcc: <0.5 0 0.25> <1.0 or 0.0 0.5 0.75>
+     -> diag: <0.5 0 0.25>
+ *
+ * - type 2:
+ *   1st fcc: <0 0 0> <0.5 0.5 0> <0.5 0 0.5> <0 0.5 0.5>
+ *   2nd fcc: <0.25 0.25 0.25> <0.75 0.75 0.25> <0.75 0.25 0.75>
+ *            <0.25 0.75 0.75>
+ *   -> diag: <0.25 0.25 0.25>
+ *
+ */
 
 int main(int argc,char **argv) {
 
-       int i,j,k,l,cnt;
+       int i,j,k,l,cnt,estimated;
+       int x,y,z;
        t_3dvec basis[3];
-       t_3dvec offset;
-       t_3dvec r,n;
-       char type;
-       double zrot;
-
-       if(argc<5) {
-               printf("usage: %s type lx ly lz\n",argv[0]);
+       t_3dvec o[3];
+       t_3dvec dia;
+       t_3dvec r,h;
+       char type,fccdia;
+
+       if(argc!=6) {
+               printf("usage: %s type fcc/dia lx ly lz\n",argv[0]);
+               printf("  basis types: 0, 1, 2 (see code)\n");
+               printf("  fcc/dia: 0=fcc, 1=dia\n");
                return -1;
        }
 
        type=argv[1][0];
-       type=DIAMOND;
+       fccdia=argv[2][0];
+
+       x=atoi(argv[3]);
+       y=atoi(argv[4]);
+       z=atoi(argv[5]);
 
-       zrot=0.0;
-       if(argc==6)
-               zrot=atof(argv[5]);
+       dia.x=0.25;
+       dia.y=0.25;
+       dia.z=0.25;
 
-       offset.x=0.25;
-       offset.y=0.25;
-       offset.z=0.25;
+       if(type=='1') {
+               dia.x=0.5;
+               dia.y=0.0;
+               dia.z=0.25;
+       }
 
        memset(basis,0,3*sizeof(t_3dvec));
-       basis[0].x=0.5;
-       basis[0].y=0.5;
-       basis[1].x=0.5;
-       basis[1].z=0.5;
-       basis[2].y=0.5;
-       basis[2].z=0.5;
+       memset(o,0,3*sizeof(t_3dvec));
+
+       if(type=='0') {
+               basis[0].x=0.5;
+               basis[0].y=0.5;
+               basis[1].y=0.5;
+               basis[1].z=0.5;
+               basis[2].x=0.5;
+               basis[2].z=0.5;
+       }
+
+       if(type=='1') {
+               basis[0].x=0.5;
+               basis[0].y=-0.5;
+               basis[1].x=0.5;
+               basis[1].y=0.5;
+               basis[2].z=1.0;
+               o[0].x=0.5;
+               o[0].y=0.5;
+               o[0].z=0.5;
+       }
+
+       if(type=='2') {
+               basis[0].x=1.0;
+               basis[1].y=1.0;
+               basis[2].z=1.0;
+               o[0].x=0.5;
+               o[0].y=0.5;
+               o[1].x=0.5;
+               o[1].z=0.5;
+               o[2].y=0.5;
+               o[2].z=0.5;
+       }
 
        cnt=0;
-       r.x=0.0;
-       for(i=0;i<atoi(argv[2]);i++) {
-               r.y=0.0;
-               for(j=0;j<atoi(argv[3]);j++) {
-                       r.z=0.0;
-                       for(k=0;k<atoi(argv[4]);k++) {
-
-       // first atom
-       printf(" %.2f %.2f %.2f T T T\n",r.x,r.y,r.z);
+
+       estimated=1;
+       if(fccdia=='1') estimated*=2;
+       if(type=='1') estimated*=2;
+       if(type=='2') estimated*=4;
+       estimated*=x*y*z;
+
+       // print POSCAR 'header'
+
+       printf("cubic diamond\n");
+       printf(" 5.429\n");
+
+       v3_scale(&h,&basis[0],x);
+       printf(" %.5f %.5f %.5f\n",h.x,h.y,h.z);
+       v3_scale(&h,&basis[1],y);
+       printf(" %.5f %.5f %.5f\n",h.x,h.y,h.z);
+       v3_scale(&h,&basis[2],z);
+       printf(" %.5f %.5f %.5f\n",h.x,h.y,h.z);
+
+       printf(" %d\n",estimated);
+       printf("selective dynamics\n");
+       printf("direct\n");
+
+       // now print the coordinates
+       
+       for(i=0;i<x;i++) {
+               for(j=0;j<y;j++) {
+                       for(k=0;k<z;k++) {
+
+       // first atom (all types)
+       r.x=i;
+       r.y=j;
+       r.z=k;
+       
+       printf(" %.5f %.5f %.5f T T T\n",r.x/x,r.y/y,r.z/z);
        cnt+=1;
-       // face centered atoms
-       for(l=0;l<3;l++) {
-               v3_add(&n,&r,&(basis[l]));
-               printf(" %.2f %.2f %.2f T T T\n",n.x,n.y,n.z);
+
+       // second atom (type 1)
+       if(type=='1') {
+               v3_add(&h,&r,&(o[0]));
+               printf(" %.5f %.5f %.5f T T T\n",h.x/x,h.y/y,h.z/z);
                cnt+=1;
        }
 
-                               r.z+=1.0;
+       // face centered atoms (type 2)
+       if(type=='2') {
+               for(l=0;l<3;l++) {
+                       v3_add(&h,&r,&(o[l]));
+                       printf(" %.5f %.5f %.5f T T T\n",h.x/x,h.y/y,h.z/z);
+                       cnt+=1;
+               }
+       }
+
                        }
-                       r.y+=1.0;
                }
-               r.x+=1.0;
        }
 
        // second fcc lattice
-       if(type==DIAMOND) {
+       if(fccdia=='1') {
+
+       for(i=0;i<x;i++) {
+               for(j=0;j<y;j++) {
+                       for(k=0;k<z;k++) {
 
-       r.x=0.0;
-       for(i=0;i<atoi(argv[2]);i++) {
-               r.y=0.0;
-               for(j=0;j<atoi(argv[3]);j++) {
-                       r.z=0.0;
-                       for(k=0;k<atoi(argv[4]);k++) {
+       // first atom (all types)
+       r.x=i;
+       r.y=j;
+       r.z=k;
 
-       // first atom
-       printf(" %.2f %.2f %.2f T T T\n",r.x+0.25,r.y+0.25,r.z+0.25);
+       // add the diagonal
+       v3_add(&r,&r,&dia);
+       
+       printf(" %.5f %.5f %.5f T T T\n",r.x/x,r.y/y,r.z/z);
        cnt+=1;
-       // face centered atoms
-       for(l=0;l<3;l++) {
-               v3_add(&n,&r,&(basis[l]));
-               printf(" %.2f %.2f %.2f T T T\n",n.x+0.25,n.y+0.25,n.z+0.25);
+
+       // second atom (type 1)
+       if(type=='1') {
+               v3_add(&h,&r,&(o[0]));
+               printf(" %.5f %.5f %.5f T T T\n",h.x/x,h.y/y,h.z/z);
                cnt+=1;
        }
 
-                               r.z+=1.0;
+       // face centered atoms (type 2)
+       if(type=='2') {
+               for(l=0;l<3;l++) {
+                       v3_add(&h,&r,&(o[l]));
+                       printf(" %.5f %.5f %.5f T T T\n",h.x/x,h.y/y,h.z/z);
+                       cnt+=1;
+               }
+       }
+
                        }
-                       r.y+=1.0;
                }
-               r.x+=1.0;
        }
 
+
        }
 
 
-       printf("\ntotal: %d ions\n\n",cnt);
+       if(estimated!=cnt) {
+               printf("\nWARNING! counted %d ions and estimated %d ions\n\n",
+                      cnt,estimated);
+       }
 
        return 0;
-
 }
-