added report system, still messing around with virial pressure
[physik/posic.git] / sic.c
1 /*
2  * sic.c - investigation of the sic precipitation process of silicon carbide
3  *
4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
5  *
6  */
7
8 #include <math.h>
9  
10 #include "moldyn.h"
11
12 #include "posic.h"
13
14 int hook(void *moldyn,void *hook_params) {
15
16         t_moldyn *md;
17
18         md=moldyn;
19
20         /* switch to direct scaling in first hook */    
21         if(md->schedule.count==0)
22                 set_pt_scale(md,0,0,T_SCALE_BERENDSEN,100.0);
23         /* switch off temp scaling in second hook */
24         if(md->schedule.count==1)
25                 set_pt_scale(md,0,0,0,0);
26                 
27         //set_temperature(md,md->t_ref-100.0);
28
29         return 0;
30 }
31
32 int main(int argc,char **argv) {
33
34         /* check argv */
35         if(argc!=3) {
36                 printf("[sic] usage: %s <logdir> <temperatur>\n",argv[0]);
37                 return -1;
38         }
39
40         /* main moldyn structure */
41         t_moldyn md;
42
43         /* potential parameters */
44         t_lj_params lj;
45         t_ho_params ho;
46         t_tersoff_mult_params tp;
47
48         /* testing location & velocity vector */
49         t_3dvec r,v;
50
51         /* initialize moldyn */
52         moldyn_init(&md,argc,argv);
53
54         /* choose integration algorithm */
55         set_int_alg(&md,MOLDYN_INTEGRATE_VERLET);
56
57         /* choose potential */
58         //set_potential1b(&md,tersoff_mult_1bp,&tp);
59         //set_potential2b(&md,tersoff_mult_2bp,&tp);
60         //set_potential2b_post(&md,tersoff_mult_post_2bp,&tp);
61         //set_potential3b(&md,tersoff_mult_3bp,&tp);
62         set_potential2b(&md,lennard_jones,&lj);
63         //set_potential2b(&md,harmonic_oscillator,&ho);
64
65         /* cutoff radius */
66         //set_cutoff(&md,TM_S_SI);
67         set_cutoff(&md,2*LC_SI*0.5*sqrt(1.5));
68
69         /*
70          * potential parameters
71          */
72
73         /* lennard jones */
74         lj.sigma6=LJ_SIGMA_SI*LJ_SIGMA_SI*LJ_SIGMA_SI;
75         lj.sigma6*=lj.sigma6;
76         lj.sigma12=lj.sigma6*lj.sigma6;
77         lj.epsilon4=4.0*LJ_EPSILON_SI;
78         lj.uc=lj.epsilon4*(lj.sigma12/pow(md.cutoff,12.0)-lj.sigma6/pow(md.cutoff,6));
79
80         /* harmonic oscillator */
81         ho.equilibrium_distance=0.25*sqrt(3.0)*LC_SI;
82         ho.spring_constant=.1;
83
84         /*
85          * tersoff mult potential parameters for SiC
86          */
87         tp.S[0]=TM_S_SI;
88         tp.R[0]=TM_R_SI;
89         tp.A[0]=TM_A_SI;
90         tp.B[0]=TM_B_SI;
91         tp.lambda[0]=TM_LAMBDA_SI;
92         tp.mu[0]=TM_MU_SI;
93         tp.beta[0]=TM_BETA_SI;
94         tp.n[0]=TM_N_SI;
95         tp.c[0]=TM_C_SI;
96         tp.d[0]=TM_D_SI;
97         tp.h[0]=TM_H_SI;
98
99         tp.S[1]=TM_S_C;
100         tp.R[1]=TM_R_C;
101         tp.A[1]=TM_A_C;
102         tp.B[1]=TM_B_C;
103         tp.lambda[1]=TM_LAMBDA_C;
104         tp.mu[1]=TM_MU_C;
105         tp.beta[1]=TM_BETA_C;
106         tp.n[1]=TM_N_C;
107         tp.c[1]=TM_C_C;
108         tp.d[1]=TM_D_C;
109         tp.h[1]=TM_H_C;
110
111         tp.chi=TM_CHI_SIC;
112
113         tersoff_mult_complete_params(&tp);
114
115         /* set (initial) dimensions of simulation volume */
116         set_dim(&md,8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),8*LC_SI*0.5*sqrt(1.5),TRUE);
117
118         /* set periodic boundary conditions in all directions */
119         set_pbc(&md,TRUE,TRUE,TRUE);
120
121         /* create the lattice / place atoms */
122         create_lattice(&md,FCC,LC_SI*0.5*sqrt(1.5),SI,M_SI,
123         //create_lattice(&md,DIAMOND,LC_SI,SI,M_SI,
124         //               ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
125                        ATOM_ATTR_2BP|ATOM_ATTR_HB,
126                        0,8,8,8);
127         moldyn_bc_check(&md);
128
129         /* testing configuration */
130         //r.x=0.28*sqrt(3)*LC_SI/2;     v.x=0;
131         //r.x=1.75*LC_SI;       v.x=-0.01;
132         //r.y=0;                v.y=0;
133         //r.z=0;                v.z=0;
134         //add_atom(&md,SI,M_SI,0,
135         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
136         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
137         //           &r,&v);
138         //r.x=-r.x;     v.x=-v.x;
139         //r.y=0;                v.y=0;
140         //r.z=0;                v.z=0;
141         //add_atom(&md,SI,M_SI,0,
142         //           ATOM_ATTR_1BP|ATOM_ATTR_2BP|ATOM_ATTR_3BP|ATOM_ATTR_HB,
143         //           ATOM_ATTR_2BP|ATOM_ATTR_HB,
144         //           &r,&v);
145
146         /* setting a nearest neighbour distance for the moldyn checks */
147         set_nn_dist(&md,0.25*sqrt(3.0)*LC_SI); /* diamond ! */
148
149         /* set temperature */
150         //set_temperature(&md,273.0+1410.0);
151         //set_temperature(&md,273.0+450.0);
152         //set_temperature(&md,273.0);
153         //set_temperature(&md,1.0);
154         //set_temperature(&md,0.0);
155         set_temperature(&md,atof(argv[2])+273.0);
156
157         /* set pressure */
158         set_pressure(&md,ATM);
159
160         /* set p/t scaling */
161         //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,
162         //                 T_SCALE_BERENDSEN,100.0);
163         //set_pt_scale(&md,0,0,T_SCALE_DIRECT,1.0);
164         //set_pt_scale(&md,P_SCALE_BERENDSEN,100.0,0,0);
165         
166         /* initial thermal fluctuations of particles (in equilibrium) */
167         thermal_init(&md,TRUE);
168
169         /* create the simulation schedule */
170         moldyn_add_schedule(&md,1001,1.0);
171         //moldyn_add_schedule(&md,501,1.0);
172         //moldyn_add_schedule(&md,501,1.0);
173
174         /* schedule hook function */
175         //moldyn_set_schedule_hook(&md,&hook,NULL);
176
177         /* activate logging */
178         moldyn_set_log_dir(&md,argv[1]);
179         moldyn_set_report(&md,"Frank Zirkelbach","Test 1");
180         moldyn_set_log(&md,LOG_TOTAL_ENERGY,1);
181         moldyn_set_log(&md,VISUAL_STEP,10);
182         moldyn_set_log(&md,CREATE_REPORT,0);
183
184         /*
185          * let's do the actual md algorithm now
186          *
187          * integration of newtons equations
188          */
189         moldyn_integrate(&md);
190
191         /* close */
192         moldyn_shutdown(&md);
193         
194         return 0;
195 }
196