fixed bond analyze, introduced more comfortable set potential function,

[physik/posic.git] / moldyn.c

diff --git a/moldyn.c b/moldyn.c
index 84cf700..def5079 100644 (file)
--- a/moldyn.c
+++ b/moldyn.c
@@ -20,6 +20,17 @@
 #include "moldyn.h"
 #include "report/report.h"
 
+/* potential includes */
+#include "potentials/harmonic_oscillator.h"
+#include "potentials/lennard_jones.h"
+#include "potentials/albe.h"
+#ifdef TERSOFF_ORIG
+#include "potentials/tersoff_orig.h"
+#else
+#include "potentials/tersoff.h"
+#endif
+
+
 /*
  * global variables, pse and atom colors (only needed here)
  */ 
@@ -228,58 +239,29 @@ int set_pbc(t_moldyn *moldyn,u8 x,u8 y,u8 z) {
        return 0;
 }
 
-int set_potential1b(t_moldyn *moldyn,pf_func1b func) {
-
-       moldyn->func1b=func;
-
-       return 0;
-}
-
-int set_potential2b(t_moldyn *moldyn,pf_func2b func) {
-
-       moldyn->func2b=func;
-
-       return 0;
-}
-
-int set_potential3b_j1(t_moldyn *moldyn,pf_func2b func) {
-
-       moldyn->func3b_j1=func;
-
-       return 0;
-}
-
-int set_potential3b_j2(t_moldyn *moldyn,pf_func2b func) {
-
-       moldyn->func3b_j2=func;
-
-       return 0;
-}
-
-int set_potential3b_j3(t_moldyn *moldyn,pf_func2b func) {
-
-       moldyn->func3b_j3=func;
-
-       return 0;
-}
-
-int set_potential3b_k1(t_moldyn *moldyn,pf_func3b func) {
-
-       moldyn->func3b_k1=func;
-
-       return 0;
-}
-
-int set_potential3b_k2(t_moldyn *moldyn,pf_func3b func) {
-
-       moldyn->func3b_k2=func;
-
-       return 0;
-}
-
-int set_potential_params(t_moldyn *moldyn,void *params) {
+int set_potential(t_moldyn *moldyn,u8 type) {
 
-       moldyn->pot_params=params;
+       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;
+                       // missing: check 2b bond func
+                       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->check_2b_bond=albe_mult_check_2b_bond;
+                       break;
+               default:
+                       printf("[moldyn] set potential: unknown type %02x\n",
+                              type);
+                       return -1;
+       }
 
        return 0;
 }
@@ -2508,8 +2490,8 @@ int calculate_pair_correlation_process(t_moldyn *moldyn,t_atom *itom,
        if(bc) check_per_bound(moldyn,&dist);
        d=v3_absolute_square(&dist);
 
-       /* ignore if greater or equal 2 times cutoff */
-       if(d>=4.0*moldyn->cutoff_square)
+       /* ignore if greater cutoff */
+       if(d>moldyn->cutoff_square)
                return 0;
 
        /* fill the slots */
@@ -2548,7 +2530,7 @@ int calculate_pair_correlation(t_moldyn *moldyn,double dr,void *ptr) {
        int i;
 
        pcc.dr=dr;
-       pcc.o1=2.0*moldyn->cutoff/dr;
+       pcc.o1=moldyn->cutoff/dr;
        pcc.o2=2*pcc.o1;
 
        if(pcc.o1*dr<=moldyn->cutoff)
@@ -2612,9 +2594,15 @@ int bond_analyze_process(t_moldyn *moldyn,t_atom *itom,t_atom *jtom,
        d=v3_absolute_square(&dist);
 
        /* ignore if greater or equal cutoff */
-       if(d>=moldyn->cutoff_square)
+       if(d>moldyn->cutoff_square)
                return 0;
 
+       /* check for potential bond */
+       if(moldyn->check_2b_bond(moldyn,itom,jtom,bc)==FALSE)
+               return 0;
+
+       d=sqrt(d);
+
        /* now count this bonding ... */
        ba=data;
 
@@ -2641,6 +2629,7 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
        // by now: # bonds of type 'a-4b' and 'b-4a' / # bonds total
 
        int qcnt;
+       int ccnt,cset;
        t_ba ba;
        int i;
        t_atom *atom;
@@ -2660,8 +2649,9 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
        memset(ba.bcnt,0,moldyn->count*sizeof(int));
 
        ba.tcnt=0;
-
        qcnt=0;
+       ccnt=0;
+       cset=0;
 
        atom=moldyn->atom;
 
@@ -2673,17 +2663,25 @@ int bond_analyze(t_moldyn *moldyn,double *quality) {
                                qcnt+=4;
                }
                else {
-                       if((ba.acnt[i]==4)&(ba.bcnt[i]==0))
+                       if((ba.acnt[i]==4)&(ba.bcnt[i]==0)) {
                                qcnt+=4;
+                               ccnt+=1;
+                       }
+                       cset+=1;
                }
        }
 
-printf("%d %d\n",qcnt,ba.tcnt);
-       if(quality)
-               *quality=1.0*qcnt/ba.tcnt;
-       else
-               printf("[moldyn] bond analyze: quality = %f\n",
-                      1.0*qcnt/ba.tcnt);
+       printf("[moldyn] bond analyze: c_cnt=%d | set=%d\n",ccnt,cset);
+       printf("[moldyn] bond analyze: q_cnt=%d | tot=%d\n",qcnt,ba.tcnt);
+
+       if(quality) {
+               quality[0]=1.0*ccnt/cset;
+               quality[1]=1.0*qcnt/ba.tcnt;
+       }
+       else {
+               printf("[moldyn] bond analyze: c_bnd_q=%f\n",1.0*qcnt/ba.tcnt);
+               printf("[moldyn] bond analyze:   tot_q=%f\n",1.0*qcnt/ba.tcnt);
+       }
 
        return 0;
 }