adapted to play around with new albe implementation

diff --git a/moldyn.c b/moldyn.c
index dd8a527..bf7b9f0 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -222,6 +222,7 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        moldyn->func3b_k2=tersoff_mult_3bp_k2;
                        moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
+               /*
                case MOLDYN_POTENTIAL_AM:
                        moldyn->func3b_j1=albe_mult_3bp_j1;
                        moldyn->func3b_k1=albe_mult_3bp_k1;
@@ -229,6 +230,17 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        moldyn->func3b_k2=albe_mult_3bp_k2;
                        moldyn->check_2b_bond=albe_mult_check_2b_bond;
                        break;
+               */
+               case MOLDYN_POTENTIAL_AM:
+                       moldyn->func1b=albe_mult_i0;
+                       moldyn->func2b=albe_mult_i0_j0;
+                       moldyn->func3b_0=albe_mult_i0_j0_k0;
+                       moldyn->func3b_1=albe_mult_i0_j1;
+                       moldyn->func3b_j1=albe_mult_i0_j2;
+                       moldyn->func3b_k2=albe_mult_i0_j2_k0;
+                       moldyn->func3b_j2=albe_mult_i0_j3;
+                       moldyn->check_2b_bond=albe_mult_check_2b_bond;
+                       break;
                case MOLDYN_POTENTIAL_HO:
                        moldyn->func2b=harmonic_oscillator;
                        moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
@@ -1901,6 +1913,9 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        if(jtom==&(itom[i]))
                                                continue;
 
+                                       /* reset 3bp run */
+                                       moldyn->run3bp=1;
+
                                        if((jtom->attr&ATOM_ATTR_2BP)&
                                           (itom[i].attr&ATOM_ATTR_2BP)) {
                                                moldyn->func2b(moldyn,
@@ -1908,6 +1923,78 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                               jtom,
                                                               bc_ij);
                                        }
+
+// REWRITE ONCE WORKING!!!
+                                       /* 3 body potential/force */
+
+                                       /* in j loop, 3bp run can be skipped */
+                                       if(!(moldyn->run3bp))
+                                               continue;
+
+                                       if(!(itom[i].attr&ATOM_ATTR_3BP))
+                                               continue;
+
+                                       if(!(jtom->attr&ATOM_ATTR_3BP))
+                                               continue;
+
+                                       if(moldyn->func3b_0==NULL)
+                                               continue;
+
+                               for(k=0;k<27;k++) {
+
+                                       bc_ik=(k<dnlc)?0:1;
+#ifdef STATIC_LISTS
+                                       q=0;
+
+                                       while(neighbour_i[j][q]!=0) {
+
+                                               ktom=&(atom[neighbour_i[k][q]]);
+                                               q++;
+#else
+                                       that=&(neighbour_i2[k]);
+                                       list_reset_f(that);
+                                       
+                                       if(that->start==NULL)
+                                               continue;
+
+                                       do {
+                                               ktom=that->current->data;
+#endif
+
+                                               if(!(ktom->attr&ATOM_ATTR_3BP))
+                                                       continue;
+
+                                               if(ktom==jtom)
+                                                       continue;
+
+                                               if(ktom==&(itom[i]))
+                                                       continue;
+
+                                               moldyn->func3b_0(moldyn,
+                                                                &(itom[i]),
+                                                                jtom,
+                                                                ktom,
+                                                                bc_ik|bc_ij);
+#ifdef STATIC_LISTS
+                                       }
+#else
+                                       } while(list_next_f(that)!=\
+                                               L_NO_NEXT_ELEMENT);
+#endif
+
+                               }
+
+
+                               if(moldyn->func3b_1)
+                                       moldyn->func3b_1(moldyn,
+                                                        &(itom[i]),
+                                                        jtom,
+                                                        bc_ij);
+
+// UNTIL HERE + change function names!
+
+
+
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 

diff --git a/moldyn.h b/moldyn.h
index ca40103..6345f3f 100644 (file)
--- a/moldyn.h
+++ b/moldyn.h
@@ -106,6 +106,9 @@ typedef struct s_moldyn {
        /* potential force function and parameter pointers */
        int (*func1b)(struct s_moldyn *moldyn,t_atom *ai);
        int (*func2b)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
+       int (*func3b_0)(struct s_moldyn *moldyn,
+                       t_atom *ai,t_atom *aj,t_atom *ak,u8 bck);
+       int (*func3b_1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        int (*func3b_j1)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        int (*func3b_j2)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);
        int (*func3b_j3)(struct s_moldyn *moldyn,t_atom *ai,t_atom *aj,u8 bc);