new albe potential + new force calc routine (old potentials need to be
[physik/posic.git] / moldyn.c
index bf7b9f0..1aff5a8 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -215,11 +215,11 @@ int set_potential(t_moldyn *moldyn,u8 type) {
 
        switch(type) {
                case MOLDYN_POTENTIAL_TM:
-                       moldyn->func1b=tersoff_mult_1bp;
-                       moldyn->func3b_j1=tersoff_mult_3bp_j1;
-                       moldyn->func3b_k1=tersoff_mult_3bp_k1;
-                       moldyn->func3b_j2=tersoff_mult_3bp_j2;
-                       moldyn->func3b_k2=tersoff_mult_3bp_k2;
+                       moldyn->func_i0=tersoff_mult_1bp;
+                       moldyn->func_j1=tersoff_mult_3bp_j1;
+                       moldyn->func_j1_k0=tersoff_mult_3bp_k1;
+                       moldyn->func_j1c=tersoff_mult_3bp_j2;
+                       moldyn->func_j1_k1=tersoff_mult_3bp_k2;
                        moldyn->check_2b_bond=tersoff_mult_check_2b_bond;
                        break;
                /*
@@ -232,21 +232,21 @@ int set_potential(t_moldyn *moldyn,u8 type) {
                        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->func_i0=albe_mult_i0;
+                       moldyn->func_j0=albe_mult_i0_j0;
+                       moldyn->func_j0_k0=albe_mult_i0_j0_k0;
+                       moldyn->func_j0e=albe_mult_i0_j1;
+                       moldyn->func_j1=albe_mult_i0_j2;
+                       moldyn->func_j1_k0=albe_mult_i0_j2_k0;
+                       moldyn->func_j1c=albe_mult_i0_j3;
                        moldyn->check_2b_bond=albe_mult_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_HO:
-                       moldyn->func2b=harmonic_oscillator;
+                       moldyn->func_j0=harmonic_oscillator;
                        moldyn->check_2b_bond=harmonic_oscillator_check_2b_bond;
                        break;
                case MOLDYN_POTENTIAL_LJ:
-                       moldyn->func2b=lennard_jones;
+                       moldyn->func_j0=lennard_jones;
                        moldyn->check_2b_bond=lennard_jones_check_2b_bond;
                        break;
                default:
@@ -1860,8 +1860,8 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                /* single particle potential/force */
                if(itom[i].attr&ATOM_ATTR_1BP)
-                       if(moldyn->func1b)
-                               moldyn->func1b(moldyn,&(itom[i]));
+                       if(moldyn->func_i0)
+                               moldyn->func_i0(moldyn,&(itom[i]));
 
                if(!(itom[i].attr&(ATOM_ATTR_2BP|ATOM_ATTR_3BP)))
                        continue;
@@ -1876,8 +1876,14 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                dnlc=lc->dnlc;
 
+#ifndef STATIC_LISTS
+               /* check for later 3 body interaction */
+               if(itom[i].attr&ATOM_ATTR_3BP)
+                       memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
+#endif
+
                /* first loop over atoms j */
-               if(moldyn->func2b) {
+               if(moldyn->func_j0) {
                        for(j=0;j<27;j++) {
 
                                bc_ij=(j<dnlc)?0:1;
@@ -1888,18 +1894,6 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        jtom=&(atom[neighbour_i[j][p]]);
                                        p++;
-
-                                       if(jtom==&(itom[i]))
-                                               continue;
-
-                                       if((jtom->attr&ATOM_ATTR_2BP)&
-                                          (itom[i].attr&ATOM_ATTR_2BP)) {
-                                               moldyn->func2b(moldyn,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
-                                       }
-                               }
 #else
                                this=&(neighbour_i[j]);
                                list_reset_f(this);
@@ -1909,6 +1903,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                do {
                                        jtom=this->current->data;
+#endif
 
                                        if(jtom==&(itom[i]))
                                                continue;
@@ -1918,13 +1913,12 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                        if((jtom->attr&ATOM_ATTR_2BP)&
                                           (itom[i].attr&ATOM_ATTR_2BP)) {
-                                               moldyn->func2b(moldyn,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
+                                               moldyn->func_j0(moldyn,
+                                                               &(itom[i]),
+                                                               jtom,
+                                                               bc_ij);
                                        }
 
-// REWRITE ONCE WORKING!!!
                                        /* 3 body potential/force */
 
                                        /* in j loop, 3bp run can be skipped */
@@ -1937,9 +1931,10 @@ int potential_force_calc(t_moldyn *moldyn) {
                                        if(!(jtom->attr&ATOM_ATTR_3BP))
                                                continue;
 
-                                       if(moldyn->func3b_0==NULL)
+                                       if(moldyn->func_j0_k0==NULL)
                                                continue;
 
+                               /* first loop over atoms k in first j loop */
                                for(k=0;k<27;k++) {
 
                                        bc_ik=(k<dnlc)?0:1;
@@ -1964,17 +1959,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_0(moldyn,
-                                                                &(itom[i]),
-                                                                jtom,
-                                                                ktom,
-                                                                bc_ik|bc_ij);
+                                               moldyn->func_j0_k0(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 #ifdef STATIC_LISTS
                                        }
 #else
@@ -1984,35 +1979,27 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                }
 
-
-                               if(moldyn->func3b_1)
-                                       moldyn->func3b_1(moldyn,
+                               /* finish of first j loop after first k loop */
+                               if(moldyn->func_j0e)
+                                       moldyn->func_j0e(moldyn,
                                                         &(itom[i]),
                                                         jtom,
                                                         bc_ij);
 
-// UNTIL HERE + change function names!
-
-
-
+#ifdef STATIC_LISTS
+                               }
+#else
                                } while(list_next_f(this)!=L_NO_NEXT_ELEMENT);
 #endif
 
                        }
                }
 
-               /* 3 body potential/force */
+               /* continued 3 body potential/force */
 
                if(!(itom[i].attr&ATOM_ATTR_3BP))
                        continue;
 
-               /* copy the neighbour lists */
-#ifdef STATIC_LISTS
-               /* no copy needed for static lists */
-#else
-               memcpy(neighbour_i2,neighbour_i,27*sizeof(t_list));
-#endif
-
                /* second loop over atoms j */
                for(j=0;j<27;j++) {
 
@@ -2045,18 +2032,18 @@ int potential_force_calc(t_moldyn *moldyn) {
                                /* reset 3bp run */
                                moldyn->run3bp=1;
 
-                               if(moldyn->func3b_j1)
-                                       moldyn->func3b_j1(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,
-                                                         bc_ij);
+                               if(moldyn->func_j1)
+                                       moldyn->func_j1(moldyn,
+                                                       &(itom[i]),
+                                                       jtom,
+                                                       bc_ij);
 
-                               /* in first j loop, 3bp run can be skipped */
+                               /* in j loop, 3bp run can be skipped */
                                if(!(moldyn->run3bp))
                                        continue;
                        
-                               /* first loop over atoms k */
-                               if(moldyn->func3b_k1) {
+                               /* first loop over atoms k in second j loop */
+                               if(moldyn->func_j1_k0) {
 
                                for(k=0;k<27;k++) {
 
@@ -2082,17 +2069,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_k1(moldyn,
-                                                                 &(itom[i]),
-                                                                 jtom,
-                                                                 ktom,
-                                                                 bc_ik|bc_ij);
+                                               moldyn->func_j1_k0(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 #ifdef STATIC_LISTS
                                        }
 #else
@@ -2104,14 +2091,15 @@ int potential_force_calc(t_moldyn *moldyn) {
 
                                }
 
-                               if(moldyn->func3b_j2)
-                                       moldyn->func3b_j2(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,
-                                                         bc_ij);
+                               /* continued j loop after first k loop */
+                               if(moldyn->func_j1c)
+                                       moldyn->func_j1c(moldyn,
+                                                        &(itom[i]),
+                                                        jtom,
+                                                        bc_ij);
 
                                /* second loop over atoms k */
-                               if(moldyn->func3b_k2) {
+                               if(moldyn->func_j1_k1) {
 
                                for(k=0;k<27;k++) {
 
@@ -2137,17 +2125,17 @@ int potential_force_calc(t_moldyn *moldyn) {
                                                if(!(ktom->attr&ATOM_ATTR_3BP))
                                                        continue;
 
-                                               if(ktom==jtom)
-                                                       continue;
+                                               //if(ktom==jtom)
+                                               //      continue;
 
                                                if(ktom==&(itom[i]))
                                                        continue;
 
-                                               moldyn->func3b_k2(moldyn,
-                                                                 &(itom[i]),
-                                                                 jtom,
-                                                                 ktom,
-                                                                 bc_ik|bc_ij);
+                                               moldyn->func_j1_k1(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
 
 #ifdef STATIC_LISTS
                                        }
@@ -2160,11 +2148,11 @@ int potential_force_calc(t_moldyn *moldyn) {
                                
                                }
 
-                               /* 2bp post function */
-                               if(moldyn->func3b_j3) {
-                                       moldyn->func3b_j3(moldyn,
-                                                         &(itom[i]),
-                                                         jtom,bc_ij);
+                               /* finish of j loop after second k loop */
+                               if(moldyn->func_j1e) {
+                                       moldyn->func_j1e(moldyn,
+                                                        &(itom[i]),
+                                                        jtom,bc_ij);
                                }
 #ifdef STATIC_LISTS
                        }
@@ -2173,7 +2161,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
                
                }
-               
+
 #ifdef DEBUG
        //printf("\n\n");
 #endif