X-Git-Url: https://hackdaworld.org/gitweb/?p=physik%2Fposic.git;a=blobdiff_plain;f=moldyn.h;h=0670f3d19239ada190d6c7d65fd99317392af295;hp=3016bc491cf1049942aafaeb8dd7e2c0b8641960;hb=HEAD;hpb=dc70c570abec4596355df26ff19756658e33e762 diff --git a/moldyn.h b/moldyn.h index 3016bc4..0670f3d 100644 --- a/moldyn.h +++ b/moldyn.h @@ -12,157 +12,539 @@ #include "random/random.h" #include "list/list.h" -/* datatypes */ +/* + * + * datatypes + * + */ +/* general */ +typedef unsigned char u8; + +/* virial */ +typedef struct s_virial { + double xx; /* | xx xy xz | */ + double yy; /* V = | yx yy yz | */ + double zz; /* | zx zy zz | */ + double xy; /* */ + double xz; /* with: xy=yx, xz=zx, yz=zy */ + double yz; /* */ +} t_virial; + +/* the atom of the md simulation */ typedef struct s_atom { - t_3dvec r; /* positions */ - t_3dvec v; /* velocities */ - t_3dvec f; /* forces */ - int element; /* number of element in pse */ - double mass; /* atom mass */ + 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 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 { - int listfd; - int nx,ny,nz; - int cells; - double x,y,z; - t_list *subcell; + int nx,ny,nz; /* amount of cells in x, y and z direction */ + 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 -typedef struct s_moldyn { - /* atoms, amount, dimensions */ +/* moldyn schedule structure */ +typedef struct s_moldyn_schedule { int count; - t_atom *atom; - t_3dvec dim; - /* potential, force & parameters */ - int (*potential_force_function)(struct s_moldyn *moldyn); + int total_sched; + int *runs; + double *tau; + int (*hook)(void *moldyn,void *hook_params); + 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 (*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; - /* cut off radius */ - double cutoff; - double cutoff_square; - /* linked list / cell method */ - t_linkcell lc; - /* temperature */ - double t; - /* integration of newtons equations */ + unsigned char run3bp; + + double cutoff; /* cutoff radius */ + double cutoff_square; /* square of the cutoff radius */ + double nnd; /* nearest neighbour distance (optional) */ + + 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 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 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.) */ + unsigned char pt_scale; /* type of p and t scaling */ + double t_tc; /* t berendsen control time constant */ + double p_tc; /* p berendsen control time constant */ + + /* simulation schedule */ + t_moldyn_schedule schedule; + int current; /* current position in schedule */ + + /* integration function pointer */ int (*integrate)(struct s_moldyn *moldyn); - int time_steps; - double tau; - double tau_square; + int time_steps; /* amount of iterations */ + double tau; /* delta t */ + double time; /* absolute time */ + double tau_square; /* delta t squared */ + int total_steps; /* total steps */ + /* energy */ - double energy; - /* logging & visualization */ - t_visual vis; - unsigned char lvstat; - unsigned int ewrite; - char efb[64]; - int efd; - unsigned int mwrite; - char mfb[64]; - int mfd; - unsigned int swrite; - char sfb[64]; - int sfd; - unsigned int dwrite; - char dfb[64]; - int dfd; - unsigned int vwrite; - char vfb[64]; - void *visual; - /* moldyn general status */ - unsigned char status; - /* random */ - t_random random; + 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 */ + unsigned int ewrite; /* how often to log energy */ + int efd; /* fd for energy log */ + unsigned int mwrite; /* how often to log momentum */ + int mfd; /* fd for momentum log */ + unsigned int pwrite; /* how often to log pressure */ + 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 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 */ + char rauthor[64]; /* report author */ + int epfd; /* energy gnuplot script file descriptor */ + int ppfd; /* pressure gnuplot script file descriptor */ + int tpfd; /* temperature gnuplot script file descriptor */ + + u8 status; /* general moldyn properties */ + + 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_ho_params { - double spring_constant; - double equilibrium_distance; -} t_ho_params; +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; -typedef struct s_lj_params { - double sigma6; - double sigma12; - double epsilon4; -} t_lj_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 + * */ -/* general defines */ +#define MOLDYN_STAT_PBX 0x01 /* periodic boudaries in x */ +#define MOLDYN_STAT_PBY 0x02 /* y */ +#define MOLDYN_STAT_PBZ 0x04 /* and z direction */ + +#define MOLDYN_PSCALE 0x08 /* size controlled by piston */ -#define MOLDYN_LVSTAT_TOTAL_E 0x01 -#define MOLDYN_LVSTAT_TOTAL_M 0x02 -#define MOLDYN_LVSTAT_SAVE 0x04 -#define MOLDYN_LVSTAT_DUMP 0x08 -#define MOLDYN_LVSTAT_VISUAL 0x10 -#define MOLDYN_LVSTAT_INITIALIZED 0x80 +#define MOLDYN_1BP 0x10 /* care about single */ +#define MOLDYN_2BP 0x20 /* 2 body */ +#define MOLDYN_3BP 0x40 /* and 3 body particle pots */ -#define MOLDYN_STAT_POTENTIAL 0x01 -#define MOLDYN_STAT_FORCE 0x02 +#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 + * + * - length unit: 1 A (1 A = 1e-10 m) + * - time unit: 1 fs (1 fs = 1e-15 s) + * - mass unit: 1 amu (1 amu = 1.6605388628e-27 kg ) + * + * fyi: in the following 1 N = (amu*A)/(fs*fs) + * + */ + +#define METER 1e10 /* A */ +#define SECOND 1e15 /* fs */ +#define AMU 1.6605388628e-27 /* kg */ +#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.0e-15 -#define MOLDYN_CUTOFF 10.0e-9 +#define MOLDYN_TAU 1.0 +#define MOLDYN_CUTOFF 10.0 #define MOLDYN_RUNS 1000000 #define MOLDYN_INTEGRATE_VERLET 0x00 -#define MOLDYN_INTEGRATE_DEFAULT MOLDYN_INTEGRATE_VERLET +#define MOLDYN_INTEGRATE_DEFAULT MOLDYN_INTEGRATE_VERLET #define MOLDYN_POTENTIAL_HO 0x00 #define MOLDYN_POTENTIAL_LJ 0x01 -#define MOLDYN_POTENTIAL_DEFAULT MOLDYN_POTENTIAL_LJ +#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 LOG_VOLUME 0x10 +#define SAVE_STEP 0x20 +#define VISUAL_STEP 0x40 +#define CREATE_REPORT 0x80 -/* phsical values */ +#define TRUE 1 +#define FALSE 0 -#define K_BOLTZMANN 1.3807e-27 /* Nm/K */ -#define AMU 1.660540e-27 /* kg */ +#define VERBOSE 1 +#define QUIET 0 -#define FCC 0x01 -#define DIAMOND 0x02 +#define SCALE_UP 'u' +#define SCALE_DOWN 'd' +#define SCALE_DIRECT 'D' + +/* + * usefull constants + */ + +#define ONE_THIRD (1.0/3.0) + +/* + * element specific defines + */ #define C 0x06 -#define M_C (12.011*AMU) +#define LC_C 3.567 /* A */ +#define M_C 12.011 /* amu */ #define SI 0x0e -#define LC_SI 0.543105e-9 /* m */ -#define M_SI (28.085*AMU) /* kg */ -#define LJ_SIGMA_SI ((0.20*sqrt(3.0)*LC_SI)/1.122462) /* m */ -#define LJ_EPSILON_SI (2.1678*1.60e-19) /* Nm */ +#define LC_SI 5.43105 /* A */ +#define M_SI 28.08553 /* amu */ + +#define LC_3C_SIC 4.3596 /* A */ -/* function prototypes */ +/* + * lattice types + */ + +#define CUBIC 0x01 +#define FCC 0x02 +#define DIAMOND 0x04 +#define ZINCBLENDE 0x08 +#define NONE 0x80 + +/* + * more includes + */ + +//#include "pse.h" + +/* + * + * function prototypes + * + */ -int moldyn_usage(char **argv); -int moldyn_parse_argv(t_moldyn *moldyn,int argc,char **argv); -int moldyn_log_init(t_moldyn *moldyn); int moldyn_init(t_moldyn *moldyn,int argc,char **argv); int moldyn_shutdown(t_moldyn *moldyn); -int create_lattice(unsigned char type,int element,double mass,double lc, - int a,int b,int c,t_atom **atom); -int destroy_lattice(t_atom *atom); -int thermal_init(t_moldyn *moldyn); -int scale_velocity(t_moldyn *moldyn); -double get_e_kin(t_atom *atom,int count); -double get_e_pot(t_moldyn *moldyn); +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_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, + 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); +int scale_dim(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z); +int scale_atoms(t_moldyn *moldyn,u8 dir,double scale,u8 x,u8 y,u8 z); + +double e_kin_calc(t_moldyn *moldyn); double get_total_energy(t_moldyn *moldyn); -t_3dvec get_total_p(t_atom *atom,int count); +t_3dvec get_total_p(t_moldyn *moldyn); -double estimate_time_step(t_moldyn *moldyn,double nn_dist,double t); +double estimate_time_step(t_moldyn *moldyn,double nn_dist); -int link_cell_init(t_moldyn *moldyn); +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 *); + +int moldyn_add_schedule(t_moldyn *moldyn,int runs,double tau); +int moldyn_set_schedule_hook(t_moldyn *moldyn,set_hook hook,void *hook_params); + int moldyn_integrate(t_moldyn *moldyn); int velocity_verlet(t_moldyn *moldyn); -int harmonic_oscillator(t_moldyn *moldyn); -int lennard_jones(t_moldyn *moldyn); +int potential_force_calc(t_moldyn *moldyn); +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 +