]> www.hackdaworld.org Git - physik/posic.git/blob - init/init.c
07b4dc6fd061e120f2af3b749b7ca2ff7462e97c
[physik/posic.git] / init / init.c
1 /*
2  * init.c - create initial conditions
3  *
4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
5  *
6  */
7
8 #include "../math/math.h"
9 #include "../moldyn.h"
10
11 /* fcc lattice init */
12 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
13
14         int count;
15         int i,j;
16         t_3dvec o,r,n;
17         t_3dvec basis[3];
18         double help[3];
19         double x,y,z;
20
21         x=a*lc;
22         y=b*lc;
23         z=c*lc;
24
25         if(origin) v3_copy(&o,origin);
26         else v3_zero(&o);
27
28         /* construct the basis */
29         for(i=0;i<3;i++) {
30                 for(j=0;j<3;j++) {
31                         if(i!=j) help[j]=0.5*lc;
32                         else help[j]=.0;
33                 }
34                 v3_set(&basis[i],help);
35         }
36
37         v3_zero(&r);
38         v3_add(&r,&r,&o);
39         count=0;
40         
41         /* fill up the room */
42         while(r.x<x) {
43                 r.y=.0;
44                 while(r.y<y) {
45                         r.z=.0;
46                         while(r.z<z) {
47                                 v3_copy(&(atom[count].r),&r);
48                                 count+=1;
49                                 for(i=0;i<3;i++) {
50                                         v3_add(&n,&r,&basis[i]);
51                                         if((n.x<x+o.x)&&(n.y<y+o.y)&&(n.z<z+o.z)) {
52                                                 v3_copy(&(atom[count].r),&n);
53                                                 count+=1;
54                                         }
55                                 }
56                                 r.z+=lc;        
57                         }
58                         r.y+=lc;
59                 }
60                 r.x+=lc;
61         }
62
63         /* coordinate transformation */
64         help[0]=x/2.0;
65         help[1]=y/2.0;
66         help[2]=z/2.0;
67         v3_set(&n,help);
68         for(i=0;i<count;i++)
69                 v3_sub(&(atom[i].r),&(atom[i].r),&n);
70                 
71         return count;
72 }
73
74 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin) {
75
76         int count;
77         t_3dvec o;
78
79         count=fcc_init(a,b,c,lc,atom,origin);
80
81         o.x=0.25*lc;
82         o.y=0.25*lc;
83         o.z=0.25*lc;
84
85         if(origin) v3_add(&o,&o,origin);
86
87         count+=fcc_init(a,b,c,lc,&atom[count],&o);
88
89         return count;
90 }