From: hackbard Date: Thu, 30 Mar 2006 13:49:53 +0000 (+0000) Subject: started potential and force stuff X-Git-Url: https://hackdaworld.org/gitweb/?a=commitdiff_plain;h=85abe46fecc79a3d7885ff6124186b2be2ca96ac;p=physik%2Fposic.git started potential and force stuff --- diff --git a/math/math.c b/math/math.c index 9cfaff6..9217c88 100644 --- a/math/math.c +++ b/math/math.c @@ -76,3 +76,20 @@ double v3_norm(t_3dvec *a) { return(sqrt(v3_absolute_square(a))); } +int v3_per_bounds(t_3dvec *a,t_3dvec *dim) { + + double x,y,z; + + x=0.5*dim->x; + y=0.5*dim->y; + z=0.5*dim->z; + + if(a->x>x) a->x-=dim->x; + else if(-a->x>x) a->x+=dim->x; + if(a->y>y) a->y-=dim->y; + else if(-a->y>y) a->y+=dim->y; + if(a->z>z) a->z-=dim->z; + else if(-a->z>z) a->z+=dim->z; + + return 0; +} diff --git a/math/math.h b/math/math.h index 53e2114..5ef7163 100644 --- a/math/math.h +++ b/math/math.h @@ -27,5 +27,6 @@ int v3_copy(t_3dvec *trg,t_3dvec *src); int v3_cmp(t_3dvec *a,t_3dvec *b); double v3_absolute_square(t_3dvec *a); double v3_norm(t_3dvec *a); +int v3_per_bounds(t_3dvec *a,t_3dvec *dim); #endif diff --git a/moldyn.c b/moldyn.c index ecba8ba..82749b2 100644 --- a/moldyn.c +++ b/moldyn.c @@ -147,6 +147,21 @@ double get_e_kin(t_atom *atom,int count) { return e; } +double get_e_pot(t_moldyn *moldyn) { + + return(moldyn->potential(moldyn)); +} + +double get_total_energy(t_moldyn *moldyn) { + + double e; + + e=get_e_kin(moldyn->atom,moldyn->count); + e+=get_e_pot(moldyn); + + return e; +} + t_3dvec get_total_p(t_atom *atom, int count) { t_3dvec p,p_total; @@ -161,7 +176,16 @@ t_3dvec get_total_p(t_atom *atom, int count) { return p_total; } -double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) { + +/* + * + * potentials & corresponding forces + * + */ + +/* lennard jones potential & force for one sort of atoms */ + +double potential_lennard_jones(t_moldyn *moldyn) { t_lj_params *params; t_atom *atom; @@ -172,7 +196,7 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) { double u; double eps,sig6,sig12; - params=ptr; + params=moldyn->pot_params; atom=moldyn->atom; count=moldyn->count; eps=params->epsilon; @@ -183,7 +207,7 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) { for(i=0;iatom; + count=moldyn->count; + params=moldyn->pot_params; + + for(i=0;idim)); + d=v3_absolute_square(&distance); + if(d<=moldyn->cutoff_square) { + h1=1.0/d; /* 1/r^2 */ + d=h1*h1; /* 1/r^4 */ + h2=d*d; /* 1/r^8 */ + h1*=d; /* 1/r^6 */ + h1*=h2; /* 1/r^14 */ + } + } + } + + return 0; +} diff --git a/moldyn.h b/moldyn.h index 30c4bf8..b43ed8e 100644 --- a/moldyn.h +++ b/moldyn.h @@ -21,12 +21,14 @@ typedef struct s_atom { double mass; /* atom mass */ } t_atom; - typedef struct s_moldyn { int count; t_atom *atom; - double (*potential)(void *ptr); - int (*force)(struct s_moldyn *moldyn,void *ptr); + double (*potential)(struct s_moldyn *moldyn); + void *pot_params; + int (*force)(struct s_moldyn *moldyn); + double cutoff_square; + t_3dvec dim; unsigned char status; } t_moldyn; @@ -68,8 +70,10 @@ int destroy_lattice(t_atom *atom); int thermal_init(t_atom *atom,t_random *random,int count,double t); int scale_velocity(t_atom *atom,int count,double t); double get_e_kin(t_atom *atom,int count); +double get_e_pot(t_moldyn *moldyn); +double get_total_energy(t_moldyn *moldyn); t_3dvec get_total_p(t_atom *atom,int count); -double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr); +double potential_lennard_jones(t_moldyn *moldyn); #endif diff --git a/posic.c b/posic.c index 4eb2d04..6d901b4 100644 --- a/posic.c +++ b/posic.c @@ -14,6 +14,8 @@ int main(int argc,char **argv) { + t_moldyn md; + t_atom *si; t_visual vis; @@ -21,10 +23,13 @@ int main(int argc,char **argv) { t_random random; int a,b,c; - double t,e; + double t,e,u; + double help; t_3dvec p; int count; + t_lj_params lj; + char fb[32]="saves/fcc_test"; /* init */ @@ -66,6 +71,30 @@ int main(int argc,char **argv) { p=get_total_p(si,count); printf("total momentum: %f\n",v3_norm(&p)); + /* check potential energy */ + md.count=count; + md.atom=si; + md.potential=potential_lennard_jones; + md.pot_params=&lj; + md.force=NULL; + md.status=0; + + lj.sigma6=3.0/16.0*LC_SI*LC_SI; + help=lj.sigma6*lj.sigma6; + lj.sigma6*=help; + lj.sigma12=lj.sigma6*lj.sigma6; + lj.epsilon=1; + + u=get_e_pot(&md); + + printf("potential energy: %f\n",u); + printf("total energy (1): %f\n",e+u); + printf("total energy (2): %f\n",get_total_energy(&md)); + + md.dim.x=a*LC_SI; + md.dim.y=b*LC_SI; + md.dim.z=c*LC_SI; + /* * let's do the actual md algorithm now *