X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.h;h=0670f3d19239ada190d6c7d65fd99317392af295;hp=5160668f5379b72a1cee78db4aa31f7ca381299c;hb=HEAD;hpb=92ef07d77a4c879527180224acea73a3f6564497 diff --git a/moldyn.h b/moldyn.h index 5160668..0670f3d 100644 --- a/moldyn.h +++ b/moldyn.h @@ -33,24 +33,37 @@ typedef struct s_virial { /* the atom of the md simulation */ typedef struct s_atom { + t_3dvec r_0; /* initial position */ t_3dvec r; /* position */ t_3dvec v; /* velocity */ t_3dvec f; /* force */ t_virial virial; /* virial */ double e; /* site energy */ + double ekin; /* kinetic energy */ int element; /* number of element in pse */ double mass; /* atom mass */ u8 brand; /* brand id */ int tag; /* atom unique id (number of atom) */ u8 attr; /* attributes */ + int pbc[3]; /* pb crossing in x, y and z direction */ } t_atom; #define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ #define ATOM_ATTR_HB 0x02 /* coupled to heat bath (velocity scaling) */ +#define ATOM_ATTR_VA 0x04 /* visualize this atom */ // TODO +#define ATOM_ATTR_VB 0x08 /* visualize the bond of this atom */ + +#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ +#define ATOM_ATTR_2BP 0x20 /* pair potential */ +#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ -#define ATOM_ATTR_1BP 0x10 /* single paricle potential */ -#define ATOM_ATTR_2BP 0x20 /* pair potential */ -#define ATOM_ATTR_3BP 0x40 /* 3 body potential */ +#define DEFAULT_ATOM_ATTR 0x74 // 1,2,3 body interaction + visualize + +/* special list structure for low mem approach */ +typedef struct s_lowmem_list { + int *head; + int *list; +} t_lowmem_list; /* cell lists */ typedef struct s_linkcell { @@ -58,11 +71,17 @@ typedef struct s_linkcell { int cells; /* total amount of cells */ double len; /* prefered cell edge length */ double x,y,z; /* the actual cell lengthes */ +#ifdef STATIC_LISTS + int **subcell; /* pointer to the cell lists */ +#elif LOWMEM_LISTS + t_lowmem_list *subcell; /* low mem approach list */ +#else t_list *subcell; /* pointer to the cell lists */ +#endif int dnlc; /* direct neighbour lists counter */ } t_linkcell; -#include "visual/visual.h" +#define MAX_ATOMS_PER_LIST 20 /* moldyn schedule structure */ typedef struct s_moldyn_schedule { @@ -74,24 +93,38 @@ typedef struct s_moldyn_schedule { void *hook_params; } t_moldyn_schedule; +/* visualization structure */ +typedef struct s_visual { + int fd; /* rasmol script file descriptor */ + char fb[128]; /* basename of the save files */ + t_3dvec dim; /* dimensions of the simulation cell */ +} t_visual; + /* moldyn main structure */ typedef struct s_moldyn { + int argc; /* number of arguments */ + char **args; /* pointer to arguments */ + int count; /* total amount of atoms */ + double mass; /* total system mass */ t_atom *atom; /* pointer to the atoms */ t_3dvec dim; /* dimensions of the simulation volume */ double volume; /* volume of sim cell (dim.x*dim.y*dim.z) */ /* potential force function and parameter pointers */ - int (*func1b)(struct s_moldyn *moldyn,t_atom *ai); - int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); - int (*func3b_k1)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); - int (*func3b_k2)(struct s_moldyn *moldyn, - t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_i0)(struct s_moldyn *moldyn,t_atom *ai); + int (*func_j0)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j0_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j0e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k0)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1c)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); + int (*func_j1_k1)(struct s_moldyn *moldyn, + t_atom *ai,t_atom *aj,t_atom *ak,u8 bck); + int (*func_j1e)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc); void *pot_params; unsigned char run3bp; @@ -101,25 +134,45 @@ typedef struct s_moldyn { t_linkcell lc; /* linked cell list interface */ + int avg_skip; /* amount of steps without average calc */ + double t_ref; /* reference temperature */ double t; /* actual temperature */ double t_sum; /* sum over all t */ - double mean_t; /* mean value of t */ - - t_virial virial; /* global virial (absolute coordinates) */ - double gp; /* pressure computed from global virial */ - double gp_sum; /* sum over all gp */ - double mean_gp; /* mean value of gp */ - - double mean_v; /* mean of virial */ + double t_avg; /* average value of t */ + + /* for sale! */ + t_virial gvir; /* global virial (absolute coordinates) */ + //double gv; + //double gv_sum; + //double gv_avg; + double sale1; + double sale2; + double sale3; + + // gp stuff exchanged by kinetic energies + //double gp; /* pressure computed from global virial */ + //double gp_sum; /* sum over all gp */ + //double gp_avg; /* average value of gp */ + double ekinx; + double ekiny; + double ekinz; + + t_virial vir; /* actual virial */ + double virial; double virial_sum; /* sum over all calculated virials */ + double virial_avg; /* average of virial */ double p_ref; /* reference pressure */ double p; /* actual pressure (computed by virial) */ + double px,py,pz; /* components of pressure */ double p_sum; /* sum over all p */ - double mean_p; /* mean value of p */ - t_3dvec tp; /* thermodynamic pressure dU/dV */ - double dv; /* dV for thermodynamic pressure calc */ + double p_avg; /* average value of p */ + + double tp; /* thermodynamic pressure dU/dV */ + double tp_sum; /* sum over dU/dV pressure */ + double tp_avg; /* average value of dU/dV pressure */ + int tp_cnt; /* how often to do thermodynamic p calc */ /* pressure and temperature control (velocity/volume scaling) */ /* (t_tc in units of tau, p_tc in units of tau * isoth. compressib.) */ @@ -139,9 +192,26 @@ typedef struct s_moldyn { double tau_square; /* delta t squared */ int total_steps; /* total steps */ + /* energy */ double energy; /* potential energy */ double ekin; /* kinetic energy */ + /* energy averages & fluctuations */ + double k_sum; /* sum of kinetic energy */ + double v_sum; /* sum of potential energy */ + double k_avg; /* average of kinetic energy */ + double v_avg; /* average of potential energy */ + double k2_sum; /* sum of kinetic energy squared */ + double v2_sum; /* sum of potential energy squared */ + double k2_avg; /* average of kinetic energy squared */ + double v2_avg; /* average of potential energy squared */ + double dk2_avg; /* mean square kinetic energy fluctuations */ + double dv2_avg; /* mean square potential energy fluctuations */ + + /* response functions */ + double c_v_nve; /* constant volume heat capacity (nve) */ + double c_v_nvt; /* constant volume heat capacity (nvt) */ + char vlsdir[128]; /* visualization/log/save directory */ t_visual vis; /* visualization interface structure */ u8 vlsprop; /* log/vis/save properties */ @@ -153,7 +223,9 @@ typedef struct s_moldyn { int pfd; /* fd for pressure log */ unsigned int twrite; /* how often to log temperature */ int tfd; /* fd for temperature log */ - unsigned int vwrite; /* how often to visualize atom information */ + unsigned int vwrite; /* how often to log volume */ + int vfd; /* fd for volume log */ + unsigned int awrite; /* how often to visualize atom information */ unsigned int swrite; /* how often to create a save file */ int rfd; /* report file descriptor */ char rtitle[64]; /* report title */ @@ -167,8 +239,66 @@ typedef struct s_moldyn { t_random random; /* random interface */ double debug; /* debugging stuff, ignore */ + + /* potential 2 body check function */ + int (*check_2b_bond)(struct s_moldyn *moldyn, + t_atom *itom,t_atom *jtom,u8 bc); } t_moldyn; +typedef struct s_pcc { + int o1; + int o2; + double dr; + double *stat; +} t_pcc; + +typedef struct s_ba { + int *acnt; + int *bcnt; + int tcnt; +} t_ba; + +typedef struct s_vb { + int fd; +} t_vb; + +typedef struct s_part_params { + u8 type; + double r; + t_3dvec p; + t_3dvec d; +} t_part_params; + +#define PART_INSIDE_R 1 +#define PART_OUTSIDE_R 2 +#define PART_INSIDE_D 3 +#define PART_OUTSIDE_D 4 + +typedef struct s_defect_params { + u8 type; + u8 stype; + double od; + double dd; + int element; + u8 brand; + u8 attr; +} t_defect_params; + +#define DEFECT_TYPE_0D 1 +#define DEFECT_TYPE_1D 2 +#define DEFECT_TYPE_2D 3 +#define DEFECT_TYPE_3D 4 + +#define DEFECT_STYPE_DB_X 1 +#define DEFECT_STYPE_DB_Y 2 +#define DEFECT_STYPE_DB_Z 3 +#define DEFECT_STYPE_DB_R 4 + +typedef struct s_offset_params { + t_3dvec o; + u8 use; +} t_offset_params; + /* * * defines @@ -185,10 +315,15 @@ typedef struct s_moldyn { #define MOLDYN_2BP 0x20 /* 2 body */ #define MOLDYN_3BP 0x40 /* and 3 body particle pots */ +#define T_SCALE_NONE 0x00 #define T_SCALE_BERENDSEN 0x01 /* berendsen t control */ #define T_SCALE_DIRECT 0x02 /* direct t control */ +#define T_SCALE_MASK 0x03 + +#define P_SCALE_NONE 0x00 #define P_SCALE_BERENDSEN 0x04 /* berendsen p control */ #define P_SCALE_DIRECT 0x08 /* direct p control */ +#define P_SCALE_MASK 0x0c /* * default values & units @@ -207,9 +342,12 @@ typedef struct s_moldyn { #define KILOGRAM (1.0/AMU) /* amu */ #define NEWTON (METER*KILOGRAM/(SECOND*SECOND)) /* A amu / fs^2 */ #define PASCAL (NEWTON/(METER*METER)) /* N / A^2 */ +#define GPA (1e9*PASCAL) /* N / A^2 */ #define BAR ((1.0e5*PASCAL)) /* N / A^2 */ #define K_BOLTZMANN (1.380650524e-23*METER*NEWTON) /* NA/K */ +#define K_B2 (K_BOLTZMANN*K_BOLTZMANN) /* (NA)^2/K^2 */ #define EV (1.6021765314e-19*METER*NEWTON) /* NA */ +#define JOULE (NEWTON*METER) /* NA */ #define MOLDYN_TEMP 273.0 #define MOLDYN_TAU 1.0 @@ -222,14 +360,17 @@ typedef struct s_moldyn { #define MOLDYN_POTENTIAL_HO 0x00 #define MOLDYN_POTENTIAL_LJ 0x01 #define MOLDYN_POTENTIAL_TM 0x02 +#define MOLDYN_POTENTIAL_AM 0x03 +#define MOLDYN_POTENTIAL_AO 0x04 #define LOG_TOTAL_ENERGY 0x01 #define LOG_TOTAL_MOMENTUM 0x02 #define LOG_PRESSURE 0x04 #define LOG_TEMPERATURE 0x08 -#define SAVE_STEP 0x10 -#define VISUAL_STEP 0x20 -#define CREATE_REPORT 0x40 +#define LOG_VOLUME 0x10 +#define SAVE_STEP 0x20 +#define VISUAL_STEP 0x40 +#define CREATE_REPORT 0x80 #define TRUE 1 #define FALSE 0 @@ -242,12 +383,15 @@ typedef struct s_moldyn { #define SCALE_DIRECT 'D' /* - * potential related phsical values / constants - * + * usefull constants */ #define ONE_THIRD (1.0/3.0) +/* + * element specific defines + */ + #define C 0x06 #define LC_C 3.567 /* A */ #define M_C 12.011 /* amu */ @@ -258,82 +402,6 @@ typedef struct s_moldyn { #define LC_3C_SIC 4.3596 /* A */ -#define LJ_SIGMA_SI ((0.25*sqrt(3.0)*LC_SI)/1.122462) /* A */ -//#define LJ_SIGMA_SI (LC_SI/1.122462) /* A */ -//#define LJ_SIGMA_SI (0.5*sqrt(2.0)*LC_SI/1.122462) /* A */ -#define LJ_EPSILON_SI (2.1678*EV) /* NA */ - -#define TM_R_SI 2.7 /* A */ -#define TM_S_SI 3.0 /* A */ -#define TM_A_SI (1830.8*EV) /* NA */ -#define TM_B_SI (471.18*EV) /* NA */ -#define TM_LAMBDA_SI 2.4799 /* 1/A */ -#define TM_MU_SI 1.7322 /* 1/A */ -#define TM_BETA_SI 1.1000e-6 -#define TM_N_SI 0.78734 -#define TM_C_SI 1.0039e5 -#define TM_D_SI 16.217 -#define TM_H_SI -0.59825 - -#define TM_R_C 1.8 /* A */ -#define TM_S_C 2.1 /* A */ -#define TM_A_C (1393.6*EV) /* NA */ -#define TM_B_C (346.7*EV) /* NA */ -#define TM_LAMBDA_C 3.4879 /* 1/A */ -#define TM_MU_C 2.2119 /* 1/A */ -#define TM_BETA_C 1.5724e-7 -#define TM_N_C 0.72751 -#define TM_C_C 3.8049e4 -#define TM_D_C 4.384 -#define TM_H_C -0.57058 - -#define TM_CHI_SIC 0.9776 - -#define TM_LC_3C_SIC (0.432e-9*METER) /* A */ - -#define ALBE_R_SI (2.82-0.14) -#define ALBE_S_SI (2.82+0.14) -#define ALBE_A_SI (3.24*EV/0.842) -#define ALBE_B_SI (1.842*3.24*EV/0.842) -#define ALBE_R0_SI 2.232 -#define ALBE_LAMBDA_SI (1.4761*sqrt(2.0*1.842)) -#define ALBE_MU_SI (1.4761*sqrt(2.0/1.842)) -#define ALBE_GAMMA_SI 0.114354 -#define ALBE_C_SI 2.00494 -#define ALBE_D_SI 0.81472 -#define ALBE_H_SI 0.259 - -#define LC_SI_ALBE 5.429 - -#define ALBE_R_C (2.00-0.15) -#define ALBE_S_C (2.00+0.15) -#define ALBE_A_C (6.00*EV/1.167) -#define ALBE_B_C (2.167*6.00*EV/1.167) -#define ALBE_R0_C 1.4276 -#define ALBE_LAMBDA_C (2.0099*sqrt(2.0*2.167)) -#define ALBE_MU_C (2.0099*sqrt(2.0/2.167)) -#define ALBE_GAMMA_C 0.11233 -#define ALBE_C_C 181.910 -#define ALBE_D_C 6.28433 -#define ALBE_H_C 0.5556 - -#define LC_C_ALBE 3.566 - -#define ALBE_R_SIC (2.40-0.20) -#define ALBE_S_SIC (2.40+0.10) -#define ALBE_A_SIC (4.36*EV/0.847) -#define ALBE_B_SIC (1.847*4.36*EV/0.847) -#define ALBE_R0_SIC 1.79 -#define ALBE_LAMBDA_SIC (1.6991*sqrt(2.0*1.847)) -#define ALBE_MU_SIC (1.6991*sqrt(2.0/1.847)) -#define ALBE_GAMMA_SIC 0.011877 -#define ALBE_C_SIC 273987 -#define ALBE_D_SIC 180.314 -#define ALBE_H_SIC 0.68 - -#define LC_SIC_ALBE 4.359 - - /* * lattice types */ @@ -341,7 +409,14 @@ typedef struct s_moldyn { #define CUBIC 0x01 #define FCC 0x02 #define DIAMOND 0x04 +#define ZINCBLENDE 0x08 +#define NONE 0x80 + +/* + * more includes + */ +//#include "pse.h" /* * @@ -349,10 +424,6 @@ typedef struct s_moldyn { * */ -typedef int (*pf_func1b)(t_moldyn *,t_atom *); -typedef int (*pf_func2b)(t_moldyn *,t_atom *,t_atom *,u8); -typedef int (*pf_func3b)(t_moldyn *,t_atom *,t_atom *,t_atom *,u8); - int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); @@ -360,38 +431,46 @@ int set_int_alg(t_moldyn *moldyn,u8 algo); int set_cutoff(t_moldyn *moldyn,double cutoff); int set_temperature(t_moldyn *moldyn,double t_ref); int set_pressure(t_moldyn *moldyn,double p_ref); +int set_p_scale(t_moldyn *moldyn,u8 ptype,double ptc); +int set_t_scale(t_moldyn *moldyn,u8 ttype,double ttc); int set_pt_scale(t_moldyn *moldyn,u8 ptype,double ptc,u8 ttype,double ttc); int set_dim(t_moldyn *moldyn,double x,double y,double z,u8 visualize); int set_nn_dist(t_moldyn *moldyn,double dist); int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z); -int set_potential1b(t_moldyn *moldyn,pf_func1b func); -int set_potential2b(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func); -int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func); -int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func); -int set_potential_params(t_moldyn *moldyn,void *params); +int set_potential(t_moldyn *moldyn,u8 type); + +int set_avg_skip(t_moldyn *moldyn,int skip); int moldyn_set_log_dir(t_moldyn *moldyn,char *dir); int moldyn_set_report(t_moldyn *moldyn,char *author,char *title); int moldyn_set_log(t_moldyn *moldyn,u8 type,int timer); 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); -int cubic_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); -int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin); -int add_atom(t_moldyn *moldyn,int element,double mass,u8 brand,u8 attr, +int create_lattice(t_moldyn *moldyn,u8 type,double lc,int element, + u8 attr,u8 brand,int a,int b,int c,t_3dvec *origin, + t_part_params *p_params,t_defect_params *d_params, + t_offset_params *o_params); +int add_atom(t_moldyn *moldyn,int element,u8 brand,u8 attr, t_3dvec *r,t_3dvec *v); +int del_atom(t_moldyn *moldyn,int tag); +int cubic_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, + t_part_params *p_params,t_defect_params *d_params); +int fcc_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, + t_part_params *p_params,t_defect_params *d_params); +int diamond_init(int a,int b,int c,double lc,t_atom *atom,t_3dvec *origin, + t_part_params *p_params,t_defect_params *d_params); int destroy_atoms(t_moldyn *moldyn); int thermal_init(t_moldyn *moldyn,u8 equi_init); +double total_mass_calc(t_moldyn *moldyn); double temperature_calc(t_moldyn *moldyn); double get_temperature(t_moldyn *moldyn); int scale_velocity(t_moldyn *moldyn,u8 equi_init); +double virial_sum(t_moldyn *moldyn); double pressure_calc(t_moldyn *moldyn); +int average_reset(t_moldyn *moldyn); +int average_and_fluctuation_calc(t_moldyn *moldyn); +int get_heat_capacity(t_moldyn *moldyn); double thermodynamic_pressure_calc(t_moldyn *moldyn); double get_pressure(t_moldyn *moldyn); int scale_volume(t_moldyn *moldyn); @@ -406,7 +485,13 @@ double estimate_time_step(t_moldyn *moldyn,double nn_dist); int link_cell_init(t_moldyn *moldyn,u8 vol); int link_cell_update(t_moldyn *moldyn); +#ifdef STATIC_LISTS +int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int **cell); +#elif LOWMEM_LISTS +int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,int *cell); +#else int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell); +#endif int link_cell_shutdown(t_moldyn *moldyn); typedef int (*set_hook)(void *,void *); @@ -422,9 +507,44 @@ int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d); //inline int virial_calc(t_atom *a,t_3dvec *f,t_3dvec *d) // __attribute__((always_inline)); int check_per_bound(t_moldyn *moldyn,t_3dvec *a); +int check_per_bound_and_care_for_pbc(t_moldyn *moldyn,t_atom *a); //inline int check_per_bound(t_moldyn *moldyn,t_3dvec *a) // __attribute__((always_inline)); int moldyn_bc_check(t_moldyn *moldyn); +int moldyn_read_save_file(t_moldyn *moldyn,char *file); +int moldyn_free_save_file(t_moldyn *moldyn); +int moldyn_load(t_moldyn *moldyn); +int process_2b_bonds(t_moldyn *moldyn,void *data, + int (*process)(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc)); +int process_neighbours(t_moldyn *moldyn,void *data,t_atom *atom, + int (*process)(t_moldyn *moldyn,t_atom *atom,t_atom *natom, + void *data,u8 bc)); + +int get_line(int fd,char *line,int max); + +int pair_correlation_init(t_moldyn *moldyn,double dr); +int calculate_diffusion_coefficient(t_moldyn *moldyn,double *dc); +int calculate_msd(t_moldyn *moldyn,double *msd); +int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom, + t_atom *jtom,void *data,u8 bc); +int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr); +int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc); +int bond_analyze(t_moldyn *moldyn,double *quality); + +int visual_init(t_moldyn *moldyn,char *filebase); +int visual_bonds_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom, + void *data,u8 bc); +#ifdef VISUAL_THREAD +void *visual_atoms(void *ptr); +#else +int visual_atoms(t_moldyn *moldyn); #endif + +int fpu_set_rtd(void); + +#endif +