From: hackbard Date: Wed, 22 Nov 2006 15:23:42 +0000 (+0000) Subject: nearly finished tersoff 2 body part and new api X-Git-Url: https://hackdaworld.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=020789bbba5881e2abf58c362d0ce733696de3b3;p=physik%2Fposic.git nearly finished tersoff 2 body part and new api --- diff --git a/moldyn.c b/moldyn.c index d74b391..9f88211 100644 --- a/moldyn.c +++ b/moldyn.c @@ -232,7 +232,7 @@ int moldyn_shutdown(t_moldyn *moldyn) { return 0; } -int create_lattice(unsigned char type,int element,double mass,double lc, +int create_lattice(u8 type,int element,double mass,double lc, int a,int b,int c,t_atom **atom) { int count; @@ -487,7 +487,7 @@ int link_cell_neighbour_index(t_moldyn *moldyn,int i,int j,int k,t_list *cell) { int ci,cj,ck; int nx,ny,nz; int x,y,z; - unsigned char bx,by,bz; + u8 bx,by,bz; lc=&(moldyn->lc); nx=lc->nx; @@ -702,7 +702,7 @@ int potential_force_calc(t_moldyn *moldyn) { t_list neighbour[27]; t_list *this; double u; - unsigned char bc,bc3; + u8 bc,bc3; int countn,dnlc; count=moldyn->count; @@ -804,304 +804,174 @@ int potential_force_calc(t_moldyn *moldyn) { return 0; } +/* + * periodic boundayr checking + */ + +int check_per_bound(t_moldyn *moldyn,t_3dvec *a) { + + double x,y,z; + + x=0.5*dim->x; + y=0.5*dim->y; + z=0.5*dim->z; + + if(moldyn->MOLDYN_ATTR_PBX) + if(a->x>=x) a->x-=dim->x; + else if(-a->x>x) a->x+=dim->x; + if(moldyn->MOLDYN_ATTR_PBY) + if(a->y>=y) a->y-=dim->y; + else if(-a->y>y) a->y+=dim->y; + if(moldyn->MOLDYN_ATTR_PBZ) + if(a->z>=z) a->z-=dim->z; + else if(-a->z>z) a->z+=dim->z; + + return 0; +} + + /* * example potentials */ /* harmonic oscillator potential and force */ -int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,unsigned char bc)) { +int harmonic_oscillator(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc)) { t_ho_params *params; - t_atom *atom,*btom; - t_linkcell *lc; - t_list *this,neighbour[27]; - int i,j,c; - int count; t_3dvec force,distance; - double d,u; + double d; double sc,equi_dist; - int ni,nj,nk; - params=moldyn->pot_params; - atom=moldyn->atom; - lc=&(moldyn->lc); + params=moldyn->pot2b_params; sc=params->spring_constant; equi_dist=params->equilibrium_distance; - count=moldyn->count; - - /* reset energy counter */ - u=0.0; - for(i=0;idim.x/2))/lc->x; - nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y; - nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z; - c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); - - /* - * processing cell of atom i - * => no need to check for empty list (1 element at minimum) - */ - this=&(neighbour[0]); - list_reset(this); - do { - btom=this->current->data; - if(btom==&(atom[i])) - continue; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - d=v3_norm(&distance); - if(d<=moldyn->cutoff) { - u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); - v3_scale(&force,&distance, - -sc*(1.0-(equi_dist/d))); - v3_add(&(atom[i].f),&(atom[i].f),&force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - - /* - * direct neighbour cells - * => no boundary condition check necessary - */ - for(j=1;jstart!=NULL) { - - do { - btom=this->current->data; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - d=v3_norm(&distance); - if(d<=moldyn->cutoff) { - u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); - v3_scale(&force,&distance, - -sc*(1.0-(equi_dist/d))); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - - } - } - - /* - * indirect neighbour cells - * => check boundary conditions - */ - for(j=c;j<27;j++) { - this=&(neighbour[j]); - list_reset(this); /* check boundary conditions */ - if(this->start!=NULL) { - - do { - btom=this->current->data; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - v3_per_bound(&distance,&(moldyn->dim)); - d=v3_norm(&distance); - if(d<=moldyn->cutoff) { - u+=(0.5*sc*(d-equi_dist)*(d-equi_dist)); - v3_scale(&force,&distance, - -sc*(1.0-(equi_dist/d))); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - - } - } + v3_sub(&distance,&(ai->r),&(aj->r); + + v3_per_bound(&distance,&(moldyn->dim)); + if(bc) check_per_bound(moldyn,&distance); + d=v3_norm(&distance); + if(d<=moldyn->cutoff) { + /* energy is 1/2 (d-d0)^2, but we will add this twice ... */ + moldyn->energy+=(0.25*sc*(d-equi_dist)*(d-equi_dist)); + v3_scale(&force,&distance,-sc*(1.0-(equi_dist/d))); + v3_add(&(ai->f),&(ai->f),&force); } - moldyn->energy=0.5*u; - return 0; } /* lennard jones potential & force for one sort of atoms */ -int lennard_jones(t_moldyn *moldyn) { +int lennard_jones(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { t_lj_params *params; - t_atom *atom,*btom; - t_linkcell *lc; - t_list *this,neighbour[27]; - int i,j,c; - int count; t_3dvec force,distance; double d,h1,h2,u; double eps,sig6,sig12; - int ni,nj,nk; params=moldyn->pot_params; - atom=moldyn->atom; - lc=&(moldyn->lc); - count=moldyn->count; eps=params->epsilon4; sig6=params->sigma6; sig12=params->sigma12; - /* reset energy counter */ - u=0.0; + v3_sub(&distance,&(ai->r),&(aj->r)); + if(bc) check_per_bound(moldyn,&distance); + d=v3_absolute_square(&distance); /* 1/r^2 */ + if(d<=moldyn->cutoff_square) { + d=1.0/d; /* 1/r^2 */ + h2=d*d; /* 1/r^4 */ + h2*=d; /* 1/r^6 */ + h1=h2*h2; /* 1/r^12 */ + /* energy is eps*..., but we will add this twice ... */ + moldyn->energy+=0.5*eps*(sig12*h1-sig6*h2); + h2*=d; /* 1/r^8 */ + h1*=d; /* 1/r^14 */ + h2*=6*sig6; + h1*=12*sig12; + d=+h1-h2; + d*=eps; + v3_scale(&force,&distance,d); + v3_add(&(ai->f),&(aj->f),&force); + } - for(i=0;idim.x/2))/lc->x; - nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y; - nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z; - c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); - - /* processing cell of atom i */ - this=&(neighbour[0]); - list_reset(this); /* list has 1 element at minimum */ - do { - btom=this->current->data; - if(btom==&(atom[i])) - continue; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - d=v3_absolute_square(&distance); /* 1/r^2 */ - if(d<=moldyn->cutoff_square) { - d=1.0/d; /* 1/r^2 */ - h2=d*d; /* 1/r^4 */ - h2*=d; /* 1/r^6 */ - h1=h2*h2; /* 1/r^12 */ - u+=eps*(sig12*h1-sig6*h2); - h2*=d; /* 1/r^8 */ - h1*=d; /* 1/r^14 */ - h2*=6*sig6; - h1*=12*sig12; - d=+h1-h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(atom[i].f),&(atom[i].f),&force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); +/* + * tersoff potential & force for 2 sorts of atoms + */ - /* neighbours not doing boundary condition overflow */ - for(j=1;jstart!=NULL) { +/* tersoff 2 body part */ - do { - btom=this->current->data; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - d=v3_absolute_square(&distance); /* r^2 */ - if(d<=moldyn->cutoff_square) { - d=1.0/d; /* 1/r^2 */ - h2=d*d; /* 1/r^4 */ - h2*=d; /* 1/r^6 */ - h1=h2*h2; /* 1/r^12 */ - u+=eps*(sig12*h1-sig6*h2); - h2*=d; /* 1/r^8 */ - h1*=d; /* 1/r^14 */ - h2*=6*sig6; - h1*=12*sig12; - d=+h1-h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); - - } - } +int tersoff_2bp(t_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc) { - /* neighbours due to boundary conditions */ - for(j=c;j<27;j++) { - this=&(neighbour[j]); - list_reset(this); /* check boundary conditions */ - if(this->start!=NULL) { + t_tersoff_params *params; + t_3dvec dist_ij; + double d_ij; - do { - btom=this->current->data; - v3_sub(&distance,&(atom[i].r),&(btom->r)); - v3_per_bound(&distance,&(moldyn->dim)); - d=v3_absolute_square(&distance); /* r^2 */ - if(d<=moldyn->cutoff_square) { - d=1.0/d; /* 1/r^2 */ - h2=d*d; /* 1/r^4 */ - h2*=d; /* 1/r^6 */ - h1=h2*h2; /* 1/r^12 */ - u+=eps*(sig12*h1-sig6*h2); - h2*=d; /* 1/r^8 */ - h1*=d; /* 1/r^14 */ - h2*=6*sig6; - h1*=12*sig12; - d=+h1-h2; - d*=eps; - v3_scale(&force,&distance,d); - v3_add(&(atom[i].f),&(atom[i].f), - &force); - } - } while(list_next(this)!=L_NO_NEXT_ELEMENT); + params=moldyn->pot_params; + + /* + * we need: f_c, df_c, f_r, df_r + * + * therefore we need: R, S + */ - } + v3_sub(&dist_ij,&(ai->r),&(aj->r)); + + if(bc) check_per_bound(moldyn,&dist_ij); + + if(ai->bnum==aj->bnum) { + S=params->S[ai->bnum]; + R=params->R[ai->bnum]; + } + else { + S=params->Smixed; + R=params->Rmixed; + } + + d_ij=v3_norm(&dist_ij); + + if(d_ij<=S) { + f_r=A*exp(-lamda*d_ij); + df_r=-lambda*f_r/d_ij; + if(d_ijf),&(ai->f),&force); } + moldyn->energy+=(f_r*f_c); } - moldyn->energy=0.5*u; - return 0; } -/* tersoff potential & force for 2 sorts of atoms */ +/* tersoff 3 body part */ -int tersoff(t_moldyn *moldyn) { +int tersoff(t_moldyn *moldyn,t_atom *ai,t_atom *aj,t_atom *ak,u8 bc,u8 bck) { t_tersoff_params *params; - t_atom *atom,*btom,*ktom; - t_linkcell *lc; - t_list *this,*thisk,neighbour[27],neighbourk[27]; - int i,j,k,c,ck; - int count; - double u; - int ni,nj,nk; - int ki,kj,kk; - + t_3dvec dist_ij; + double d_ij; params=moldyn->pot_params; - atom=moldyn->atom; - lc=&(moldyn->lc); - count=moldyn->count; - /* reset energy counter */ - u=0.0; - - for(i=0;idim.x/2))/lc->x; - nj=(atom[i].r.y+(moldyn->dim.y/2))/lc->y; - nk=(atom[i].r.z+(moldyn->dim.z/2))/lc->z; - c=link_cell_neighbour_index(moldyn,ni,nj,nk,neighbour); - - /* - * processing cell of atom i - * => no need to check for empty list (1 element at minimum) - */ - this=&(neighbour[0]); - list_reset(this); - do { - btom=this->current->data; - if(btom==&(atom[i])) - continue; - - /* 2 body stuff */ - - /* we need: f_c, df_c, f_r, df_r */ + /* 2 body part of the tersoff potential */ - v3_sub(&dist_ij,btom,&(atom[i])); - d_ij=v3_norm(&dist_ij); - if(d_ij<=S) { + v3_sub(&dist_ij,&(ai->r),&(aj->r)); + if(bc) check_per_bound(moldyn,&dist_ij); + d_ij=v3_norm(&dist_ij); + if(d_ij<=S) { /* determine the tersoff parameters */ if(atom[i].element!=btom->element) { diff --git a/moldyn.h b/moldyn.h index 66f774e..d98b428 100644 --- a/moldyn.h +++ b/moldyn.h @@ -25,7 +25,8 @@ typedef struct s_atom { t_3dvec f; /* force */ int element; /* number of element in pse */ double mass; /* atom mass */ - unsigned char attr; /* attributes */ + u8 bnum; /* brand number */ + u8 attr; /* attributes */ } t_atom; #define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ @@ -58,12 +59,14 @@ typedef struct s_moldyn { t_3dvec dim; /* dimensions of the simulation volume */ - /* potential force function pointer and parameters */ + /* potential force function and parameter pointers */ int (*pf_func1b)(struct s_moldyn *,t_atom *); + void *pot1b_params; int (*pf_func2b)(struct s_moldyn *,t_atom *,t_atom *); + void *pot2b_params; int (*pf_func3b)(struct s_moldyn *,t_atom *,t_atom *,t_atom *); + void *pot3b_params; //int (*potential_force_function)(struct s_moldyn *moldyn); - void *pot_params; /* parameters describing the potential */ double cutoff; /* cutoff radius */ double cutoff_square; /* square of the cutoff radius */ @@ -81,7 +84,7 @@ typedef struct s_moldyn { double energy; /* energy */ t_visual vis; /* visualization/log/save interface structure */ - unsigned char lvstat; /* log & vis properties */ + u8 lvstat; /* log & vis properties */ unsigned int ewrite; /* how often to log energy */ char efb[64]; /* energy log filename */ int efd; /* fd for energy log */ @@ -93,7 +96,7 @@ typedef struct s_moldyn { void *visual; /* pointer (hack!) */ unsigned int swrite; /* how often to create a save file */ - unsigned char status; /* general moldyn properties */ + u8 status; /* general moldyn properties */ t_random random; /* random interface */ } t_moldyn; @@ -120,6 +123,11 @@ typedef struct s_lj_params { } t_lj_params; typedef struct s_tersoff_params { + double S[2]; /* tersoff cutoff radii */ + double R[2]; /* tersoff cutoff radii */ + double Smixed /* mixed S radius */ + double Rmixed /* mixed R radius */ + double l_1,l_2; double m_1,m_2; double a_1,a_2; @@ -171,7 +179,7 @@ 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 create_lattice(u8 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);