Merge branch 'leadoff'

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 33224a6..1d6e0b7 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -34,6 +34,7 @@
 #include "potentials/harmonic_oscillator.h"
 #include "potentials/lennard_jones.h"
 #include "potentials/albe.h"
+#include "potentials/albe_orig.h"
 #ifdef TERSOFF_ORIG
 #include "potentials/tersoff_orig.h"
 #else
@@ -270,19 +271,29 @@ 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_AO:
+                       moldyn->func_j1=albe_orig_mult_3bp_j1;
+                       moldyn->func_j1_k0=albe_orig_mult_3bp_k1;
+                       moldyn->func_j1c=albe_orig_mult_3bp_j2;
+                       moldyn->func_j1_k1=albe_orig_mult_3bp_k2;
+                       moldyn->check_2b_bond=albe_orig_mult_check_2b_bond;
+                       break;
                case MOLDYN_POTENTIAL_AM:
-                       moldyn->func3b_j1=albe_mult_3bp_j1;
-                       moldyn->func3b_k1=albe_mult_3bp_k1;
-                       moldyn->func3b_j2=albe_mult_3bp_j2;
-                       moldyn->func3b_k2=albe_mult_3bp_k2;
+                       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:
@@ -2487,8 +2498,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;
@@ -2503,8 +2514,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;
@@ -2536,12 +2553,69 @@ 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,
-                                                              &(itom[i]),
-                                                              jtom,
-                                                              bc_ij);
+                                               moldyn->func_j0(moldyn,
+                                                               &(itom[i]),
+                                                               jtom,
+                                                               bc_ij);
+                                       }
+
+                                       /* 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->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;
+#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->func_j0_k0(moldyn,
+                                                                  &(itom[i]),
+                                                                  jtom,
+                                                                  ktom,
+                                                                  bc_ik|bc_ij);
+#ifdef STATIC_LISTS
                                        }
 #ifdef STATIC_LISTS
                                }
@@ -2554,7 +2628,7 @@ int potential_force_calc(t_moldyn *moldyn) {
                        }
                }
 
-               /* 3 body potential/force */
+               /* continued 3 body potential/force */
 
                if(!(itom[i].attr&ATOM_ATTR_3BP))
                        continue;
@@ -2607,18 +2681,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++) {
 
@@ -2651,8 +2725,8 @@ 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;
@@ -2676,14 +2750,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++) {
 
@@ -2716,17 +2791,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
                                        }
@@ -2741,11 +2816,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
                        }
@@ -2756,7 +2831,7 @@ int potential_force_calc(t_moldyn *moldyn) {
 #endif
                
                }
-               
+
 #ifdef DEBUG
        //printf("\n\n");
 #endif