From 0d71b08d4382b23b8217342caf138cfce36226cc Mon Sep 17 00:00:00 2001 From: hackbard Date: Fri, 17 Nov 2006 10:19:50 +0000 Subject: [PATCH] began rewrite ... --- moldyn.c | 24 ++++++++++++++++++------ moldyn.h | 47 +++++++++++++++++++++++++++++++---------------- 2 files changed, 49 insertions(+), 22 deletions(-) diff --git a/moldyn.c b/moldyn.c index ad0f9a3..f92ce55 100644 --- a/moldyn.c +++ b/moldyn.c @@ -481,7 +481,7 @@ int link_cell_update(t_moldyn *moldyn) { for(i=0;icells;i++) list_destroy(&(moldyn->lc.subcell[i])); - for(count=0;countcount;count++) { + for(count=0;countcount;count++) { i=(atom[count].r.x+(moldyn->dim.x/2))/lc->x; j=(atom[count].r.y+(moldyn->dim.y/2))/lc->y; k=(atom[count].r.z+(moldyn->dim.z/2))/lc->z; @@ -995,12 +995,16 @@ int tersoff(t_moldyn *moldyn) { /* 2 body stuff */ + /* we need: f_c, df_c, f_r, df_r */ + v3_sub(&dist_ij,btom,&(atom[i])); d_ij=v3_norm(&dist_ij); if(d_ij<=S) { - S=; - R=; + /* determine the tersoff parameters */ + if(atom[i].element!=btom->element) { + S=sqrt(TERSOFF_S[e1]*TERSOFF_S[e2]); + R=R_m; A=; lambda=; B=; @@ -1010,8 +1014,10 @@ int tersoff(t_moldyn *moldyn) { betaN=; if(d_ij<=R) { - f_c=1.0; - df_c=0.0; + df_r=-lambda*A*exp(-lambda*d_ij)/d_ij; + v3_scale(&force,&dist_ij,df_r); + v3_add(&(atom[i].f),&(atom[i].f), + &force); } else { s_r=S-R; @@ -1019,7 +1025,11 @@ int tersoff(t_moldyn *moldyn) { f_c=0.5+0.5*cos(arg1); df_c=-0.5*sin(arg1)*(PI/(s_r*d_ij)); f_r=A*exp(-lambda*d_ij); - f_a=-B*exp(-mu*d_ij); + df_r=-lambda*f_r/d_ij; + scale=df_c*f_r+df_r*f_c; + v3_scale(&force,&dist_ij,scale); + v3_add(&(atom[i].f),&(atom[i].f), + &force); } } else @@ -1084,6 +1094,8 @@ int tersoff(t_moldyn *moldyn) { v3_sub(&distance_jk,ktom,btom); cos_theta=(d_ij2+d_ik*d_ik-d_jk*d_jk)/\ (2*d_ij*d_ik); + sin_theta=sqrt(1.0/\ + (cos_theta*cos_theta)); theta=arccos(cos_theta); diff --git a/moldyn.h b/moldyn.h index 49ca0b2..c6712ba 100644 --- a/moldyn.h +++ b/moldyn.h @@ -12,35 +12,50 @@ #include "random/random.h" #include "list/list.h" + /* datatypes */ + +/* 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; /* position */ + t_3dvec v; /* velocity */ + t_3dvec f; /* force */ + int element; /* number of element in pse */ + double mass; /* atom mass */ + unsigned char attr; /* attributes */ } t_atom; +#define ATOM_ATTR_FP 0x01 /* fixed position (bulk material) */ +#define ATOM_ATTR_HB 0x02 /* coupled to heat bath (velocity scaling) */ + +/* 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 */ + t_list *subcell; /* pointer to the cell lists */ + int dnlc; /* direct neighbour lists counter */ } t_linkcell; #include "visual/visual.h" +# moldyn properties structure */ + typedef struct s_moldyn { - /* atoms, amount, dimensions */ - int count; - t_atom *atom; - t_3dvec dim; + + int count; /* total amount of atoms */ + t_atom *atom; /* pointer to the atoms */ + t_3dvec dim; /* dimensions of the simulation volume */ /* potential, force & parameters */ + + /* potential force funtion created by the user */ int (*potential_force_function)(struct s_moldyn *moldyn); - void *pot_params; - /* cut off radius */ + void *pot_params; /* parameters describing the potential */ + double cutoff; double cutoff_square; /* linked list / cell method */ -- 2.20.1