]> hackdaworld.org Git - physik/posic.git/blob - moldyn.h
init in moldyn now
[physik/posic.git] / moldyn.h
1 /*
2  * moldyn.h - molecular dynamics library header file
3  *
4  * author: Frank Zirkelbach <frank.zirkelbach@physik.uni-augsburg.de>
5  *
6  */
7
8 #ifndef MOLDYN_H
9 #define MOLDYN_H
10
11 #include "math/math.h"
12 #include "random/random.h"
13 #include "list/list.h"
14
15
16 /*
17  * 
18  * datatypes
19  *
20  */
21
22 /* general */
23 typedef unsigned char u8;
24
25 /* the atom of the md simulation */
26 typedef struct s_atom {
27         t_3dvec r;              /* position */
28         t_3dvec v;              /* velocity */
29         t_3dvec f;              /* force */
30         t_3dvec virial;         /* virial */
31         int element;            /* number of element in pse */
32         double mass;            /* atom mass */
33         u8 bnum;                /* brand number */
34         u8 attr;                /* attributes */
35 } t_atom;
36
37 #define ATOM_ATTR_FP    0x01    /* fixed position (bulk material) */
38 #define ATOM_ATTR_HB    0x02    /* coupled to heat bath (velocity scaling) */
39
40 #define ATOM_ATTR_1BP           0x10    /* single paricle potential */
41 #define ATOM_ATTR_2BP           0x20    /* pair potential */
42 #define ATOM_ATTR_3BP           0x40    /* 3 body potential */ 
43
44 /* cell lists */
45 typedef struct s_linkcell {
46         int nx,ny,nz;           /* amount of cells in x, y and z direction */
47         int cells;              /* total amount of cells */
48         double len;             /* prefered cell edge length */
49         double x,y,z;           /* the actual cell lengthes */
50         t_list *subcell;        /* pointer to the cell lists */
51         int dnlc;               /* direct neighbour lists counter */
52 } t_linkcell;
53
54 #include "visual/visual.h"
55
56 /* moldyn schedule structure */
57 typedef struct s_moldyn_schedule {
58         int content_count;
59         int *runs;
60         double *tau;
61         int (*hook)(void *moldyn,void *hook);
62         void *hook_params;
63 } t_moldyn_schedule;
64
65 /* moldyn main structure */
66 typedef struct s_moldyn {
67         int count;              /* total amount of atoms */
68         t_atom *atom;           /* pointer to the atoms */
69
70         t_3dvec dim;            /* dimensions of the simulation volume */
71         double volume;          /* volume of sim cell (dim.x*dim.y*dim.z) */
72
73         /* potential force function and parameter pointers */
74         int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
75         void *pot1b_params;
76         int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
77         int (*func2b_post)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
78         void *pot2b_params;
79         int (*func3b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,
80                       u8 bck);
81         void *pot3b_params;
82         //int (*potential_force_function)(struct s_moldyn *moldyn);
83
84         double cutoff;          /* cutoff radius */
85         double cutoff_square;   /* square of the cutoff radius */
86         double nnd;             /* nearest neighbour distance (optional) */
87
88         t_linkcell lc;          /* linked cell list interface */
89
90         double t_ref;           /* reference temperature */
91         double t;               /* actual temperature */
92
93         double p_ref;           /* reference pressure */
94         double p;               /* actual pressure */
95
96         /* pressure and temperature control (velocity/volume scaling) */
97         /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */
98         unsigned char pt_scale; /* type of p and t scaling */
99         double t_tc;            /* t berendsen control time constant */
100         double p_tc;            /* p berendsen control time constant */
101
102         /* simulation schedule */
103         t_moldyn_schedule schedule;
104         int current;            /* current position in schedule */
105
106         /* integration function pointer */
107         int (*integrate)(struct s_moldyn *moldyn);
108         int time_steps;         /* amount of iterations */
109         double tau;             /* delta t */
110         double time;            /* absolute time */
111         double tau_square;      /* delta t squared */
112         double elapsed;         /* total elapsed time */
113
114         double energy;          /* potential energy */
115         double ekin;            /* kinetic energy */
116
117         char vlsdir[128];       /* visualization/log/save directory */
118         t_visual vis;           /* visualization interface structure */
119         u8 vlsprop;             /* log/vis/save properties */
120         unsigned int ewrite;    /* how often to log energy */
121         int efd;                /* fd for energy log */
122         unsigned int mwrite;    /* how often to log momentum */
123         int mfd;                /* fd for momentum log */
124         unsigned int vwrite;    /* how often to visualize atom information */
125         unsigned int swrite;    /* how often to create a save file */
126
127         u8 status;              /* general moldyn properties */
128
129         t_random random;        /* random interface */
130
131         int debug;              /* debugging stuff, ignore */
132 } t_moldyn;
133
134 #define MOLDYN_STAT_PBX                 0x01    /* periodic boudaries in x */
135 #define MOLDYN_STAT_PBY                 0x02    /* y */
136 #define MOLDYN_STAT_PBZ                 0x04    /* and z direction */
137
138 #define MOLDYN_PSCALE                   0x08    /* size controlled by piston */
139
140 #define MOLDYN_1BP                      0x10    /* care about single */
141 #define MOLDYN_2BP                      0x20    /* 2 body */
142 #define MOLDYN_3BP                      0x40    /* and 3 body particle pots */
143
144 #define T_SCALE_BERENDSEN               0x01    /* berendsen t control */
145 #define T_SCALE_DIRECT                  0x02    /* direct t control */
146 #define P_SCALE_BERENDSEN               0x04    /* berendsen p control */
147 #define P_SCALE_DIRECT                  0x08    /* direct p control */
148
149
150 /*
151  *
152  * potential parameter structures
153  *
154  */
155
156 /*
157  * harmonic oscillator potential parameter structure
158  */
159
160 typedef struct s_ho_params {
161         double spring_constant;
162         double equilibrium_distance;
163 } t_ho_params;
164
165 /*
166  * lennard jones potential parameter structure
167  */
168
169 typedef struct s_lj_params {
170         double sigma6;
171         double sigma12;
172         double epsilon4;
173 } t_lj_params;
174
175 /*
176  * tersoff 
177  */
178
179 /* tersoff exchange structure to exchange 2bp and 3bp calculated values */
180 typedef struct s_tersoff_exchange {
181         double f_c,df_c;
182         double f_a,df_a;
183
184         t_3dvec dist_ij;
185         double d_ij;
186
187         double chi;
188
189         double *beta_i;
190         double *beta_j;
191         double *n_i;
192         double *n_j;
193         double *c_i;
194         double *c_j;
195         double *d_i;
196         double *d_j;
197         double *h_i;
198         double *h_j;
199
200         double ci2;
201         double cj2;
202         double di2;
203         double dj2;
204         double ci2di2;
205         double cj2dj2;
206         double betaini;
207         double betajnj;
208
209         u8 run3bp;
210         u8 run2bp_post;
211         u8 d_ij_between_rs;
212
213         double zeta_ij;
214         double zeta_ji;
215         t_3dvec dzeta_ij;
216         t_3dvec dzeta_ji;
217 } t_tersoff_exchange;
218
219 /* tersoff multi (2!) potential parameters */
220 typedef struct s_tersoff_mult_params {
221         double S[2];            /* tersoff cutoff radii */
222         double R[2];            /* tersoff cutoff radii */
223         double Smixed;          /* mixed S radius */
224         double Rmixed;          /* mixed R radius */
225         double A[2];            /* factor of tersoff attractive part */
226         double B[2];            /* factor of tersoff repulsive part */
227         double Amixed;          /* mixed A factor */
228         double Bmixed;          /* mixed B factor */
229         double lambda[2];       /* tersoff lambda */
230         double lambda_m;        /* mixed lambda */
231         double mu[2];           /* tersoff mu */
232         double mu_m;            /* mixed mu */
233
234         double chi;
235
236         double beta[2];
237         double n[2];
238         double c[2];
239         double d[2];
240         double h[2];
241
242         t_tersoff_exchange exchange;    /* exchange between 2bp and 3bp calc */
243 } t_tersoff_mult_params;
244
245
246
247 /*
248  *
249  *  defines
250  *
251  */
252
253 #define ONE_THIRD       (1.0/3.0)
254
255 /*
256  * default values
257  *
258  * - length unit: 1 A (1 A = 1e-10 m)
259  * - time unit: 1 fs (1 fs = 1e-15 s)
260  * - mass unit: 1 amu (1 amu = 1.6605388628e-27 kg )
261  *
262  * fyi: in the following 1 N = (amu*A)/(fs*fs)
263  *
264  */
265
266 #define METER                           1e10                    /* A */
267 #define SECOND                          1e15                    /* fs */
268 #define AMU                             1.6605388628e-27        /* kg */
269 #define KILOGRAM                        (1.0/AMU)               /* amu */
270 #define NEWTON  (METER*KILOGRAM/(SECOND*SECOND))        /* A amu / fs^2 */
271 #define PASCAL  (NEWTON/(METER*METER))                  /* N / A^2 */
272 #define ATM     (1.0133e5*PASCAL)                       /* N / A^2 */
273
274 #define MOLDYN_TEMP                     273.0
275 #define MOLDYN_TAU                      1.0
276 #define MOLDYN_CUTOFF                   10.0
277 #define MOLDYN_RUNS                     1000000
278
279 #define MOLDYN_INTEGRATE_VERLET         0x00
280 #define MOLDYN_INTEGRATE_DEFAULT        MOLDYN_INTEGRATE_VERLET
281
282 #define MOLDYN_POTENTIAL_HO             0x00
283 #define MOLDYN_POTENTIAL_LJ             0x01
284 #define MOLDYN_POTENTIAL_TM             0x02
285
286 #define LOG_TOTAL_ENERGY                0x01
287 #define LOG_TOTAL_MOMENTUM              0x02
288 #define SAVE_STEP                       0x04
289 #define VISUAL_STEP                     0x08
290
291 #define TRUE                            1
292 #define FALSE                           0
293
294 /*
295  *
296  * phsical values / constants
297  *
298  *
299  */
300
301 #define K_BOLTZMANN             (1.380650524e-23*METER*NEWTON)  /* NA/K */
302 #define EV                      (1.6021765314e-19*METER*NEWTON) /* NA */
303
304 #define C                       0x06
305 #define M_C                     12.011                          /* amu */
306
307 #define SI                      0x0e
308 #define LC_SI                   (0.543105e-9*METER)             /* A */
309 #define M_SI                    28.08553                        /* amu */
310 #define LJ_SIGMA_SI             ((0.25*sqrt(3.0)*LC_SI)/1.122462)       /* A */
311 #define LJ_EPSILON_SI           (2.1678*EV)                             /* NA */
312
313 #define TM_R_SI                 (2.7e-10*METER)                 /* A */
314 #define TM_S_SI                 (3.0e-10*METER)                 /* A */
315 #define TM_A_SI                 (1830.8*EV)                     /* NA */
316 #define TM_B_SI                 (471.18*EV)                     /* NA */
317 #define TM_LAMBDA_SI            (2.4799e10/METER)               /* 1/A */
318 #define TM_MU_SI                (1.7322e10/METER)               /* 1/A */
319 #define TM_BETA_SI              1.1000e-6
320 #define TM_N_SI                 0.78734
321 #define TM_C_SI                 1.0039e5
322 #define TM_D_SI                 16.217
323 #define TM_H_SI                 -0.59825
324
325 #define TM_R_C                  (1.8e-10*METER)                 /* A */
326 #define TM_S_C                  (2.1e-10*METER)                 /* A */
327 #define TM_A_C                  (1393.6*EV)                     /* NA */
328 #define TM_B_C                  (346.7*EV)                      /* NA */
329 #define TM_LAMBDA_C             (3.4879e10/METER)               /* 1/A */
330 #define TM_MU_C                 (2.2119e10/METER)               /* 1/A */
331 #define TM_BETA_C               1.5724e-7
332 #define TM_N_C                  0.72751
333 #define TM_C_C                  3.8049e4
334 #define TM_D_C                  4.384
335 #define TM_H_C                  -0.57058
336
337 #define TM_CHI_SIC              0.9776
338
339 /*
340  * lattice constants
341  */
342
343 #define FCC                     0x01
344 #define DIAMOND                 0x02
345
346
347 /*
348  *
349  * function prototypes
350  *
351  */
352
353 typedef int (*pf_func1b)(t_moldyn *,t_atom *ai);
354 typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8 bc);
355 typedef int (*pf_func2b_post)(t_moldyn *,t_atom *,t_atom *,u8 bc);
356 typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8 bc);
357
358 int moldyn_init(t_moldyn *moldyn,int argc,char **argv);
359 int moldyn_shutdown(t_moldyn *moldyn);
360
361 int set_int_alg(t_moldyn *moldyn,u8 algo);
362 int set_cutoff(t_moldyn *moldyn,double cutoff);
363 int set_temperature(t_moldyn *moldyn,double t_ref);
364 int set_pressure(t_moldyn *moldyn,double p_ref);
365 int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc);
366 int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize);
367 int set_nn_dist(t_moldyn *moldyn,double dist);
368 int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z);
369 int set_potential1b(t_moldyn *moldyn,pf_func1b func,void *params);
370 int set_potential2b(t_moldyn *moldyn,pf_func2b func,void *params);
371 int set_potential2b_post(t_moldyn *moldyn,pf_func2b_post func,void *params);
372 int set_potential3b(t_moldyn *moldyn,pf_func3b func,void *params);
373
374 int moldyn_set_log_dir(t_moldyn *moldyn,char *dir);
375 int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer);
376 int moldyn_log_shutdown(t_moldyn *moldyn);
377
378 int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element,double mass,
379                    u8 attr,u8 bnum,int a,int b,int c);
380 int add_atom(t_moldyn *moldyn,int element,double mass,u8 bnum,u8 attr,
381              t_3dvec *r,t_3dvec *v);
382 int destroy_atoms(t_moldyn *moldyn);
383
384 int thermal_init(t_moldyn *moldyn,u8 equi_init);
385 int scale_velocity(t_moldyn *moldyn,u8 equi_init);
386 int scale_volume(t_moldyn *moldyn);
387
388 double get_e_kin(t_moldyn *moldyn);
389 double get_e_pot(t_moldyn *moldyn);
390 double get_total_energy(t_moldyn *moldyn);
391 t_3dvec get_total_p(t_moldyn *moldyn);
392
393 double estimate_time_step(t_moldyn *moldyn,double nn_dist);
394
395 int link_cell_init(t_moldyn *moldyn);
396 int link_cell_update(t_moldyn *moldyn);
397 int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell);
398 int link_cell_shutdown(t_moldyn *moldyn);
399
400 int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau);
401 int moldyn_set_schedule_hook(t_moldyn *moldyn,void *hook,void *hook_params);
402
403 int moldyn_integrate(t_moldyn *moldyn);
404 int velocity_verlet(t_moldyn *moldyn);
405
406 int potential_force_calc(t_moldyn *moldyn);
407 int check_per_bound(t_moldyn *moldyn,t_3dvec *a);
408 int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
409 int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
410 int tersoff_mult_complete_params(t_tersoff_mult_params *p);
411 int tersoff_mult_1bp(t_moldyn *moldyn,t_atom *ai);
412 int tersoff_mult_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
413 int tersoff_mult_post_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
414 int tersoff_mult_3bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc);
415
416 int moldyn_bc_check(t_moldyn *moldyn);
417
418 int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
419 int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin);
420
421 #endif