started potential and force stuff

diff --git a/math/math.c b/math/math.c
index 9cfaff6..9217c88 100644 (file)
--- 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 (file)
--- 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 (file)
--- 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;i<count;i++) {
                for(j=0;j<i;j++) {
                        v3_sub(&distance,&(atom[j].r),&(atom[i].r));
-                       d=v3_absolute_square(&distance);        /* 1/r^2 */
+                       d=1.0/v3_absolute_square(&distance);    /* 1/r^2 */
                        help=d*d;                               /* 1/r^4 */
                        help*=d;                                /* 1/r^6 */
                        d=help*help;                            /* 1/r^12 */
@@ -194,4 +218,34 @@ double potential_lennard_jones_2(t_moldyn *moldyn,void *ptr) {
        return u;
 }
 
+int force_lennard_jones(t_moldyn *moldyn) {
+
+       t_lj_params *params;
+       int i,j,count;
+       t_atom *atom;
+       t_3dvec distance;
+       t_3dvec force;
+       double d,h1,h2;
+
+       atom=moldyn->atom;      
+       count=moldyn->count;
+       params=moldyn->pot_params;
+
+       for(i=0;i<count;i++) {
+               for(j=0;j<i;j++) {
+                       v3_sub(&distance,&(atom[j].r),&(atom[i].r));
+                       v3_per_bound(&distance,&(moldyn->dim));
+                       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 (file)
--- 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 (file)
--- 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
         *