started potential and force stuff

[physik/posic.git] / moldyn.c

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;
+}